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Abstract While sunspots are easily observed at the solar surface, determining 
their subsurface structure is not trivial. There are two main hypotheses for the 
subsurface structure of sunspots: the monolithic model and the cluster model. 
Local helioseismology is the only means by which we can investigate subpho- 
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tospheric structure. However, as current linear inversion techniques do not yet 
allow helioseismology to probe the internal structure with sufficient confidence to 
distinguish between the monolith and cluster models, the development of physi- 
cally realistic sunspot models are a priority for helioseismologists. This is because 
they are not only important indicators of the variety of physical effects that may 
influence helioseismic inferences in active regions, but they also enable detailed 
assessments of the validity of helioseismic interpretations through numerical 
forward modeling. In this paper, we provide a critical review of the existing 
sunspot models and an overview of numerical methods employed to model wave 
propagation through model sunspots. We then carry out an helioseismic analysis 
of the sunspot in Active Region 9787 and address the serious inconsistencies 
uncovered by iGizon et al. (2009a 2009b). We find that this sunspot is most 



probably associated with a shallow, positive wave-speed perturbation (unlike the 
traditional two-layer model) and that travel-time measurements are consistent 
with a horizontal outflow in the surrounding moat. 
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1. Introduction 



1.1. Why Sunspots are Interesting 



Sunspots are the surface manifestations of intense magnetic flux concentrations 
that have intersected with the solar surface. As such, they represent one of 
the major connections of the internal magnetic field of the Sun witli its wider 
environments and are the main sites of solar activity phenomena. They are at the 
center of many ongoing challenges in the study of the Sun, as the structure and 
evolution of sunspots, individually and collectively, are still not fully understood. 

Sunspots tend to appear at well defined latitudes, which vary with the 11- 
year solar cycle, as summarized in the so-called butterfly diagram. Any theory 
on the mechanism of the solar global dynamo has to be able to explain this 
collective behavior. Understanding dynamo processes is of utmost importance, 
as they are believed to play crucial roles in many astrophysical phenomena, and 
sunspots are the best known candidates to provide us important clues on how 
they operate. 

While the sunspots are easily observed at the surface, determining their sub- 
surface structure is not at all trivial. There are two main hypotheses for the 
structure of the subsurface magnetic configuration of the spot: the monolithic 
mod el (e.g. |Cowling|1946 1957 1976 ) and the jellyfish/cluster/spaghetti model 
{e.g. Parker|1975 1979 Spruit 1981 Zwaan 1981 ). Determining the parameters 
of these tubes, that is, typical size, field strength etc., will help reveal details 
of the operation of the solar dynamo and how magnetic field is transported up 
through the convection zone. 

Our understanding of sunspot structure has been somewhat hampered in 
the past by a lack of high-resolution observations. However, recent advances in 
adaptive optics, image selection and reconstruction at ground-based telescopes, 
and the advent of high-resolution space observations with Hinode (see Figure 
[T]), have all led to a wealth of detailed information about the fine structure of 



sunspot umbrae and penumbrae (see e.g. Scharmer et al. 2002 Thomas and 



[Weiss 2008). On the other hand, the subsurface structure of sunspots is still 
poorly understood. 

Local helioseismology {e.g. time-distance helioseismology, helioseismic holog- 
raphy, acoustic imaging, ring-diagram analysis, see [Gizon and Birch] |2005| and 
references therein) is the only means by which we can investigate subphoto- 
spheric structure. However, interpretations of data have been somewhat am- 
biguous and inconsistent for applications of local helioseismic methods in solar 
active regions. Furthermore, current linear inversion techniques do not yet allow 
helioseismology to probe the internal structure of sunspots with sufficient accu- 
racy to distinguish between monolith and cluster models. But progress has been 



made in addressing some of these inconsistencies {e.g. Braun and Birch, 2008 



Birch et al. 2009[ Gizon et al. 2009a Gizon et al.\ 2009b ) and significant ad 



vances have also been made in the simulation of helioseismic wave-propagation in 



magnetized plasmas {e.g. 


Gameron, Gizon, and DaifFallah 


2007 


Parchevsky and 


Kosovichev |20071 'Cameron, Gizon, and Duvall, 2008 ; Hanasoge 20081 Moradi, 


Hanasoge, and Gaily, 2009; Khomcnko et al. 


2009 ; Parchevsky and Kosovichev 


2009 


Shelyag et al. 


2009 


Cameron et al. 


2010 


). 
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Figure 1. A Hinode G-band image of the well-developed sunspot AR 10953, acquired at 
10:10:16UT, 2 May 2007 (the spatial scale is indicated by the horizontal bar in the upper left). 
There is a light-bridge close to the northern (upper) tip, and another is beginning to form near 
the southern edge. 



1.2. The Need for Sunspot Models in Helioseismology 

Sunspot models are important indicators of the variety of physical effects that 
may influence helioseismic inferences in active regions as well as to enable assess- 
ments of the validity of helioseismic inversions. Currently it is also important, 
in the absence of adequate nonlinear inversion techniques, to have models that 
may be close to the truth as starting points for linear inversions. The associated 
danger of course is an over-reliance on a small range of models that may limit our 
imagination of what structures may exist and which may bias the helioseismic 
inversion results. 

As illustrated in Figure[2j there are several existing classes of sunspot models - 
radiative magnetohydrodynamical (MHD) numerical simulations, semi-empirical 
models, monolithic models, cluster/spaghetti models, self-similar models, poten- 
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Figure 2. An overview diagram representing the different ciasscs of existing sunspot models. 
A number of references for each class are provided as a guide. 
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tial field models, current sheet models, pressm'e distributed models, etc. - which 
are more or less suited to local helioseismology. 

The availability of powerful computers has only recently made it possible 
to produce numerical simulations of sunspots using realistic physics (radiative 



MHD equations; see Section 4.4 1. These simulations are more than adequate for 
studying surface dynamics and fine structure, conducting experiments regarding 
the subsurface structure, and for providing very useful artificial oscillation data 
that can be used to test helioseismic inversion methods. Although this is a 
huge step forward, this approach is unfortunately not suitable for computational 
helioseismology (where detailed parametric studies need to be conducted using 
existing numerical wave propagation codes) for a number of reasons: i) The 
subsurface structure of the "realistic" simulations depends on the lower boundary 
and initial conditions, neither of which is necessarily realistic. The surface prop- 
erties obtained are more or less correct because the timescales for the evolution at 
the surface are essentially determined by the radiative loss term {i.e convection 
on the Sun is driven from the surface). The timescale for the subsurface struc- 
ture of the spot is huge, and is dependent on both the bottom the boundary 
and initial conditions. Nevertheless, work is currently underway to establish 
whether models with such artificially restrictive lower boundary conditions are 
able produce results that are compatible with helioseismic measurements, ii) 
The cumbersome computational requirements of such ab-initio simulations. As 
a guide, simulating the two-hour evolution of a pair of sunspots in |Rempel et al.\ 



(2009b) takes a number of weeks on extremely powerful supercomputers. Hi) 
Such simulations still do not address the question of the nature of the deep 
structure of sunspots, as while one may ask whether a sunspot model with the 
magnetic field clamped 10 Mm below the photosphere is, or is not, compatible 
with helioseismic measurements, probing down to these depths however will first 
require that we adequately model (and remove) the effects of the surface layers. 

Examining a number of different sunspot models is not only essential to 
test the robustness of helioseismic inferences, but it also an indispensable need 
in computational heliosiesmology where effective inference by matching wave- 
field simulations to observations relies on the availability of sunspot models. 
For example, sunspot models with very different subsurface configurations {e.g. 
cluster or monolith, connected or disconnected etc.) can be used to verify what 
helioseismic signatures would be observed for a sunspot with that particular 
subsurface structure. It is also essential to be able to continually tune the surface 
parameters of the model to actual observations in order to match the wave-field 
signatures observed. So the aim here is not to have the best sunspot model, but 
to have a group of models than can be useful for computational helioseismology - 
models which can be tuned to match the non-helioseismic observations, and then 
in turn be used in inversions to determine the subsurface structure. Of course, 
with multiple models the inversions are unlikely to be unique. Therefore, such 
parametric studies will provide us with idea of what is, and what is not reliable in 
the inversions. As well as more or less realistic models, highly simplified models 
which isolate just a few physical effects can also be useful for elucidating the sen- 
sitivity of inference methods in helioseismology to particular effects. Nonetheless, 
it must be borne in mind if using such models that other physical effects may 
have seismic signatures that are qualitatively or quantitatively similar. 
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1.3. What Basic Properties should be included in the Models? 



Apart from the desired characteristics mentioned above, models of the gross 
magnetic structure of sunspots for use in computational helioseismology should 
ideally posses a high degree of flexibility and computational efficiency to allow 
for extensive parameter studies using existing numerical simulation codes. A 
number of other essential, observationally derived, characteristics should also be 
embodied by the models. For example, an accurate prescription for the surface 
(photospheric) magnetic-field characteristics (of both the umbra and penum- 
bra, e.g. field strength, orientation, twist, return flux etc., see Section 2.2.1 and 



Section 2.2.2 1 . These should comply with observations and ideally, one should 



be able to model the sunspot field on extrapolations from observed surface 
magnetic profiles {e.g. vector magnetograms). One should also be able to choose 
the profile of thermodynamic parameters (pressure, density, temperature etc., 
see Section 2.3 and Section 4.1) in the umbra and penumbra from either spectro- 



polarimetric inversions or semi-empirical models of the solar atmosphere. Some 
important dynamical phenomena {e.g. the Evershed flow, moat flow, etc., see 



Section 2.2.4 and Section 2.2.5 ) should also be taken into account, while a realistic 



and consistent (see e.g. Section 2.4) description of the Wilson depression is also 
essential. 



1.4. The Premise of the Article 



The basic premise of this article is to satisfy two complementary goals. The first 
goal is to present a critical review of the existing physical models for the subsur- 
face structure of sunspots, in the context of local helioseismology and numerical 
simulations of wave-fields, and magnetic field-wave interaction. As discussed 
above, physical sunspot models are critically important to assess the validity of 
the helioseismic inversions. In addition, numerical simulations of the propagation 
of solar waves through model sunspots are emerging as a valid and realistic 
technique to interpret helioseismic data. The success of this approach relies on 
a very close interaction between sunspot modelers and helioseismologists. 

The second goal is to extend the helioseismic analysis undertaken by IGizon] 



et al. { 2009a 2009b,) of the sunspot in Active Region ( AR) 9787, which was the 
topic of the Third HELAS (European Hello- and Asteroseismology Network) 
Local Helioseismology Workshop, held in Berlin on 12-15 May 2009. This 
sunspot was observed during the period 20 - 28 January 2002 by the SOHO /MDI 
instrument. Serious inconsistencies between the different helioseismic methods 
were uncovered, which cannot be left unanswered. 



2. Surface Observational Constraints 

In this section we briefly review some of the main observational characteristics 
of sunspot formation and evolution. The aim here is to present a very general 
overview of some of the pertinent issues related to sunspot observations that 
need to be considered when developing a realistic sunspot model. Comprehensive 
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reviews by|Solanki|(2003) 


Thomas and Weiss ( 


2004) 


Tobias and Weiss 


2004), 


Schlichenmaier (2009 


), and the books by 


Thomas and Weiss 


(1992 


200E 


) (and 



references therein) , have excellent extended discussions on both the observational 
and theoretical aspects of sunspots and active regions. 

2.1. Sunspot Formation and Evolution 

2.1.1. Flux Emergence 

Individual bundles of magnetic flux are believed to rise from deep in the con- 
vection zone and break through the surface of the Sun. As they approach the 
surface, each flux bundle is shredded into many separate strands which, upon 
emergence, are quickly concentrated into small, intense (kG strength) magnetic 
flux bundles (or elements) by the vigorous convection occurring in the thin 
superadiabatic layer at the top of the convection zone. These small flux elements 
then accumulate at the boundaries between granules or supergranules, and some 



of them coalesce to form small pores (Keppens and Martinez Fillet, 1996 Leka 



and Skumanich, 1998). Some pores and flux elements in turn coalesce to form 
sunspots. 

A number of studies have examined the buoyant rise of magnetic flux sitting 
initially just below the surface and rising into a non-magnetized atmosphere 



[e.g. 


Magara and Longcope 




2001 




Magara and Longcope 


et al. 


2004 Cheung, Schiissler, and Moreno-Insertis 


2007 





2003 Manchester 



modeling the rise of thin flux tubes through the convection zone have also shown 
that the tubes must have a significant amount of twist in order to maintain 
their integrity and not fragment in the face of hydrodynamic forces, and indeed 
observations show that magnetic flux usually emerges at the surface already in 



a significantly twisted state {e.g. Riethmiiller et al. 2008). 



2.1.2. Sunspot Formation and Decay 

When the flux does emerge, then it is often in the form of pore structures of the 



order of a Mm or so in size ( Vrabec, 1974 Zwaan, 1978 Zwaan, 1992 Mcintosh, 



1981). Fores have continuum intensities ranging from 80% down to 20% of the 
normal photospheric intensity, with maximum magnetic-field strengths of 1500- 
2000 G. If a growing pore reaches a sufficient size (a diameter of about 3500 
km, but somet imes as much as 7000 km) or a sufficient total magnetic flux of 
order 10^° Mx (Leka and Skumanich, 1998), and if the magnetic field reaches an 
inclination from the vertical that is greater than about 7 = 35° (Martinez Fillet, 



1997), then it forms a penumbra at its periphery and becomes a fully fledged 
sunspot. The formation of a penumbra is a rapid event, occurring in less than 
20-30 minutes, and the characteristic sunspot magnetic-field configuration and 



Evershed flow are both established within this same short time period ( Leka and 



Skumanich, 1998 Yang et al., 20031. Furthermore, the fact that the largest pores 
are observed to be bigger than the smallest sunspots also provides evidence that 
the pore-sunspot transition is associated with hysteresis (Bray and Loughhead 
1964; Rucklidge, Schmidt, and Weiss 1995). 
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The time scale for the formation of a large sunspot is between a few hours 
and several days. A sunspot can span a lifetime of months, but more typically 
of weeks fSolanki, 2003). However, this life expectancy is considerably shorter 
than the magnetic diffusion time, tu = I'^/v (where / is the width of the current 
sheet and 77 = l/cj is the magnetic diffusivity) , across a solar active region where 



estimates for tu range from hundreds, to thousands of years (see Priest and 



Forbes 20001. This reduced lifetime suggests that a convective instability sets 
in that enhances the decay process, via fragmentation. Another possible process 
is the action of turbulent diffusion, owing to the nonlinear dependence of the 



diffusivity on the strength of the magnetic field (Petrovay and Moreno- Insertis, 



1997) 



An overview on sunspot decay was presented by Martinez Fillet ( 2002 1 . 
2.1.3. Surface Evolution 



Observational evidence indicates a significant change in the dynamical proper- 
ties of sunspots and active regions from "active" to "passive" evolution shortly 



after their emergence (Schiissler, 1987 Schrijver and Title, 1999 Schiissler and 
Rempel, 2005[ ). Initially, the emerging magnetic flux displays the characteristics 



of a rising, fragmented flux tube and evolves according to its internal large-scale 
dynamics {e.g. Mcintosh 1981 Strous et al. 1996). Within a few days after 
emergence, the proper motion of the sunspots with respect to the surrounding 
plasma begins to decay. Larger magnetic structures also start to fragment into 
small-scale flux concentrations, which are largely dominated by the local near- 
surface flows (granulation, supergranulation, differential rotation, meridional 
circulation). The magnetic flux is then passively advected by these velocity fields, 
gradually dispersing over a large area. This process is efficiently described by 
well-established surface fiux-transport mechanisms which model the evolution 
of the magnetic field at the solar surface {e.g. 
van Ballegooijen, Cartledge, and Priest[ [T9981 
2004f^ 



Wang, Nash, and Sheeley 



Schrijver] 2001[ Baumann et al. 



1989 



2.2. Sunspot Surface Structure 
2.2.1. Field Strength and Orientation 

A sunspot 's maximum magnetic- field strength tends to increase approximately 
linearly with sunspot diameter, from around 2000 G for the smallest, to 4000 



G or more for the largest sunspot {e.g. Ringnes and Jensen[ 
Zwaan[ [19821 |Kopp and Rabin| [1M2I [Collados et al.\ |1994| 



1960 


Brants and 


Liviui 


5ston 




2002 



periphery, becoming 700-1000 G at the edge of the visible sunspot {e.g. Mathew 



et al. 



2004). 

The field strength decreases with height within the visible outline of the spot. 
At photospheric levels in the umbra, line-of-sight field strength decreases of ~ 1 - 
3 Gkm~^ are observed (Balthasar and Schmidt, 1993 Schmidt and Balthasar, 



1994), but when averaged over a height range of 2000 km or more, the field 



strength reduces to 0.3 - 0.6 G km ^ ( Solanki, 2003 ). The strongest field within a 
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sunspot is usually associated with the darkest part of the umbra and is generally 
close to vertical. At the visible outer sunspot boundary it is inclined by 70 - 80° 
to the vertical ('Bellot Rubio et al., 2003a Mathew et al, 20041 if one applies 



inversions that assume a constant field inclination in the atmosphere. Assuming 
more than one magnetic component or gradients along the line-of-sight, high- 
resolution polarimetric studies present evidence for magnetic flux that returns 



to the surface in the outer penumbra {e.g. Westendorp Plaza et al. 2001 


Bellot 


Rubio, Balthasar, and CoUados 20041 |Borrero et al. 2004 Langhans et al. 


2005 


Borrero et al. |2006; 'Ichimoto et al. |2007| Ichimoto et al. 2009 Beck 


2008 


Jurcak and Bellot Rubio, ,2008), i.e., there is evidence that a fraction of the 



sunspot magnetic field does not extend into the chromosphere, but submerges 
beneath the photosphere in the outer penumbra. 



2.2.2. Field Twist 



As we have already mentioned, it appears that at least a small amount of twist 
is needed in a rising flux tube in order for it to resist fragmentation and preserve 
its identity as it rises through the convection zone, hence systematic twists of 
sunspot fields are relevant to the emergence and stability of sunspot magnetic 
fields. Observation also indicate that the magnetic field of regular sunspots can 
be twisted, with an azimuthal twist (j) ~ 10° -35° (Hagyard, West, and Cumings, 
1977[ Gurman and House, 1981| Lites and Skumanich, 1990| Skumanich, Lites, 
and Martinez Fillet, 19941 IWestendorp Plaza et al., 20011). 



However, the more recent observations of Mathew et al. (2003) indicate that 



for regular isolated sunspots, the global azimuthal twist of the field does not 



significantly exceed 20°. Moreover, Yun (1971) and Osherovich and Flaa (1983) 



have included the effects of an azimuthal field twist in their (self-similar) sunspot 
models and find that the introduction of a moderately twisted field, compati- 
ble with observations, contributes little to the force balance in spots and only 
slightly changes the main characteristics of their sunspot models {e.g. the mixing 
length parameter, effective temperature, Wilson depression, and the central field 
strength remain practically the same). This is not surprising, in the view of the 
fact that the measured from observations is small compared to and Br 
over most of the sunspot region. 



2.2.3. Umbral Dots 



Umbral dots, bright dot-like feature observed inside an umbra, are found in 



almost all sunspots and also in pores (Sobotka, 1997 Sobotka, 2002). They are 
observed to cover only 3-10% of the umbral area but contribute 10-20% of the 



total umbral brightness (Sobotka, Bonet, and Vazquez, 1993). Their distribution 



is not uniform: they can occur in clusters and alignments, and no large dots are 
found in dark nuclei ( [Rimmele, 1997 ). 

On average, umbral dots are 500-1000 K cooler than the photosphere outside 



a spot, but about 1000 K hotter than the coolest parts of the umbra itself (Kitai 



et al, 2007). The magnetic field in umbral dots appears to be weaker than 



that in the umbral background (Sobotka, 1997). Socas-Navarro et al. (2004) 
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found differences of several hundred gauss and deduced that the fields were more 
inclined to the vertical, by about 10° . Furthermore, a small upward velocity of 
30-50 m s~^ and 200 m has been reported in umbral dots relative to the 
umbral background by Rimmele^ (1997) and Socas-Navarro| ( ,2002) . 



Schiissler and Vogler ( 2006 1 have presented realistic numerical simulations of 



umbral magnetoconvection in the context of the monolithic model, by assuming 
an initially uniform vertical magnetic field. Their model reproduces all of the 



principal observed features of umbral dots, including their dark lanes (Rimmele, 



2008 ) . Their results provide support for the monolithic model, demonstrating 



that umbral dots can arise naturally as a consequence of magnetoconvection in 
a space-filling vertical magnetic field. The more recent numerical simulations 
of Heinemann et al. (2007) and Rempel, Schiissler, and Knolker (2009a) also 



confirm this. 



2.2.4. Ever shed Flow 



The Evershed flow is observed as an outward directed (horizontal) flow observed 
in the photospheric layers of penumbrae. The inverse Evershed flow is an inward 
directed flow in chromospheric layers. A number of well-established observational 
properties of the Evershed flow are: 

• The averaged flow velocity increases from the inner to the outer penumbra. 
From observed line asymmetries it is concluded that the flow is located 



Tritschler, 2004). 



in the deep photosphere (Maltby, 1964 Schlichenmaier, Bellot Rubio, and 



Observed velocities of the flow typically exceed 6 km s ( Rouppe van der 



Voort, 2002 Bellot Rubio et al., 2003a) in the outer penumbra, but can 
also exceed 10 km s~^ in localized patches {e.g. Bellot Rubio, Balthasar, 



and Collados 2004). 



The direction of the Evershed flow is close to horizontal. Azimuthal averages 
reveal that the flow angle varies from 70° (average flow points upwards) in 
the inner penumbra to some 100° (flow points downward) in the outer 



penumbra (Schlichenmaier and Schmidt, 2000). 



The Evershed flow is magnetized. This is obvious from Stokes V profiles 
with more than two lobes. These additional lobes are Doppler shifted {e.g. 
Schlichenmaier and Collados 2002 1 . 



Two models have been proposed to explain a number of observational prop- 



erties: the "siphon-flow" model as proposed by Montesinos and Thomas (1997), 



and the "moving-tube" model by Schlichenmaier, Jahn, and Schmidt (1998). 
Montesinos and Thomas (1997) elaborated on the idea of [Meyer and Schmidt 
( 1968 ) that the flow is driven by a gas pressure difference between the footpoints 



of a thin magnetic flux tube in magnetohydrostatic (MHS) equilibrium. On the 



other hand, Schlichenmaier, Jahn, and Schmidt ( 1998 ) developed a dynamical 



2D model of a thin magnetic flux tube that acts as a convective element in 
a superadiabatic and magnetized penumbral atmosphere. In their model, the 
convective rise of the thin flux tube to the surface initiates a local pressure 
gradient build up, leading to a gas flow along the tube. Penumbral grains are 
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then identified as the hot upfiow locations where the gas reaches the (optical) 
surface. 

Numerical simulations of radiative magnetoconvection in inclined magnetic 



fields {e.g. 


Heinemann et al. 


2007 


Scharmer, Nordlund, and Heinemann 


Rempel, Schiissler, and Knolker 2009a Rempel et al. 


2009b 


) are only beg 



2008 



to reproduce the structure of the outer penumbra, with its horizontal and return- 
ing magnetic fields and fast Evershed flows along arched channels. They have 
already succeeded in reproducing single elongated filaments with lengths of up to 
a few Mm which resemble in many ways what is observed as thin light bridges and 
penumbral filaments of the inner penumbra. In their simulations, they find that 
the progression of the filament heads toward the umbra during their formation 
phase is not caused by the inward motion of a narrow flux tube, but rather due 
to the expansion of the sheet-like upflow plumes along the filament. 

2.2.5. Moat Flow and Moving Magnetic Features 

The moat flow is an outflow that initiates immediately after the formation of 
a penumbra. Moats are typically 10 to 20 Mm wide, with the outer radius of 
the moat appearing to scale with the size of the enclosed sunspot, being about 
twice the radius of the spot itself (Brickhouse and Labonte, 19881. The moat 



flow velocity is about 0.5 to 1 km s , and can be seen by proper motions 



of granules as well as by Doppler shift measurements (Balthasar et al., 1996). 
The flow usually persists over the duration of the spot's life, while the area 
and magnetic flux of the sunspot decrease at a roughly constant rate. As the 
moat flow evolves, it pushes the magnetic flux to its periphery, leaving the moat 
largely free of magnetic field except for small magnetic features (known as moving 



km s ^ ( 


Sheeley, 1969 


Sheeley, 1972 


Vrabec, 1971 


Vrabec, 1974 Harvey and 


Harvey, 1973 Brickhouse and Labonte, 1988 Wang and Zirin, 1992 


|Yurchyshyn, 


Wang, and Goode, 2001 Sainz Dalda and Bellot Rubio, 2008aj). 



The moat flow is only present in the presence of a penumbra. Pores that 
have no penumbra also lack the moat flow. Observations of irregular sunspots 



by Vargas Dominguez et al. (2007) indicate that the moat flow exists only on 
sunspots sides where a penumbra has formed. On sunspot sides where the umbra 
and the granulation are adjacent, no moat flow is detected. Moreover, in such 
irregular configurations moat flows are only observed as radial extensions of 
penumbral filaments, but not perpendicular to the filament. 

The moat flow has also been detected through local helioseismology, using 
/-mode time-distance helioseismology. Using azimuthally-averaged MDI data. 
Gizon, Duvall, and Larsen] ( |2000 ) find an outfiow extending well beyond the 
sunspot boundary (up to 30 Mm) which reaches a peak of 1 km s~^ just outside 
the penumbra. This flow is consistent with the moat fiow. In another recent 



study, Gizon et al. ( 2009a 2009b ) used / to pi ridge-filtered time-distance travel 
times to produce linear inversion for flows around AR 9787. These inversions 
showed an azimuthally-averaged horizontal outflow in the first 4 Mm beneath 
the surface, reaching an amplitude of 230 m s~^ at a depth of 2.6 Mm and 
radial distance of some 5 Mm outside the outer spot boundary. The inversion 
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results were in line with observations of the moat flow in AR 9787 presented by 



Gizon et al. ( 2009a 2009b I , the strength and extent of which was characterized 



by measuring the observed motion of the moving magnetic features (MMFs) 
using a local correlation tracking method. Their measurements indicated a peak 
amplitude of 230 m s^^ for the moat flow, extending out to about 45 Mm from 
the spot center. 

As the moat flow is unmagnetized and has velocities that are much smaller 
than the Evershed flow, the link between the two is not obvious. One possibility 
is that the gas pressure that builds up beneath the penumbra could drive the 



moat flow (beneath the penumbra the magnetopause is largely inclined). Sainz 
Dalda and Bellot Rubio (2008b) detect bipolar MMFs within the penumbra 



which migrate outwards into and throughout the moat. This is consistent with 



an scenario proposed by Schlichenmaier ( 2002 1 : a magneto-convective overshoot 



instability in an Evershed flow channel leads to a bipolar MMF that travels 
outwards along the magnetic canopy. 



2.3. Sunspot Thermodynamics in the Photosphere 



There are a number of semi-empirical and observational models as reference, 
consisting of both one- and two-component models, for the umbra and for the 
penumbra (see Section 4.1 for a more detailed discussion). The basic assumption 
underlying almost all single component models is that it is possible to describe 
all umbrae (or at least those above a certain size) by a single thermal model. 

Thus, a very important question in this context is whether sunspot bright- 
ness (temperature) or magnetic field strength actually varies with the size of 
the sunspot. Theoretical models of sunspots based on MHS equilibrium and 



inhibition of convective heat transport {e.g. Deinzer 1965 Yun, 1970) typically 
obtain lower temperatures for stronger magnetic fields. However, observations by 



Rossbach and Schroter ( 1970 ) and Albregtsen and Maltby ( 1981 ) indicated no 



dependence of umbral intensity on umbral size for umbral diameters greater than 



about 8". On the other hand, the observations of Kopp and Rabin ( 1992 ) showed 



a nearly linear decrease in umbral brightness with umbral radius for six sunspots. 
A similar result was obtained for seven spots at visible wavelengths by |Martfnez| 
Fillet and Vazquez (1993) and Collados et al. (1994). This was also confirmed in 



a more recent study of continuum images of more than 160 sunspots taken during 



solar Cycle 23 by Mathew et al. (2007). Indeed, Norton and Oilman (2004 1 use 



this dependence to predict the peak field strength of a sunspot from its brightness 
to an accuracy of about 100 O. Hence, even though the relatively homogeneous 
umbral nuclei cannot be described by a single universal atmosphere, these models 
may nonetheless be used as boundary conditions in theoretical global sunspot 
models. 

Another related question is whether the umbra-photosphere brightness ratio 
of large sunspots varies over the solar cycle. [Albregtsen and Maltby"] ( |1978[ |1981[ ) 
find that sunspots were darkest at the beginning of sunspot Cycle 20 and that 
spots appearing later in the cycle were progressively brighter (with a nearly 
linear dependence on the phase of the cycle), up until the new Cycle 21 spots 
appeared which were again darkest. Subsequent observations showed the same 
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behavior occurred in Cycle 21 (Maltby et al, 1986). Penn and Livingston (2006) 



found similar behavior during Cycle 23. However, the results of Mathew et al.\ 
(2007) find no significant change in umbral brightness over Cycle 23. 



2.4. The Wilson Depression 

The reduced opacity in a sunspot, and the consequent depression of the t = 1 
level, arise mainly from two effects: i) the reduced temperature in the spot 
atmosphere leads to a decrease in the bound- free opacity, and ii) the radial 
force balance (including magnetic pressure and curvature forces) demands a 
lower gas pressure within the spot, further reducing the net opacity. 

A purely observational determination of the Wilson depression is complicated 
by evolutionary changes in the shape of the spot. Estimates from observations 
range from 400-1500 km for mature spots (Bray and Loughhead, 1964 Gokhale 



and Zwaan, 1972 Balthasar and Woehl, 1983). Martinez Pillet and Vazquez 



( 1993 ) assume a linear relation between magnetic pressure and temperature in a 
spot (as indicated by observations) , in which case the radial force balance yields 
a simple relationship between the net magnetic curvature force and the Wilson 
depression. They computed a Wilson depression of 400 - 800 km in the umbra of 
their observed sunspot (the range being due to different assumed values of the 
curvature forces). Solanki, Walther, and Livingston (1993) and Mathew et al. 



(2004) used this approach to determine the variation of the Wilson depression 



with radius across a sunspot. They found values of some 100 km in the penumbra 
and some 400 km in the umbra, with a fairly sharp transition at the umbra- 



penumbra boundary. In a more recent study, Watson et al. ( 2009 ) compared 



observations of sunspot longitude distribution and Monte Carlo simulations of 
sunspot appearance using different models for spot growth rate, growth time 
and depth of Wilson depression, deducing a mean depth for the umbral r = 1 
layer of 500-1500 km. 



3. The Deep Structure of Sunspots 

Though the small-scale structure seen at the surface of a sunspot is not exactly 
in equilibrium, the spot itself survives on much longer time scales than would 
be expected if it were a superficial structure. It lives for much longer than the 
time scale on which it would evolve if it were not in a stable equilibrium, which 
is the time for the Alfven speed to cross the size of the spot (on the order of an 
hour) . 

The magnetic field of a spot cannot be just a surface phenomenon, however, 
since magnetic-field lines have no ends. The extension of the spot's field lines 
above the solar surface can be observed in the chromosphere and corona, but 
they must continue below the surface as well. In contrast with a scalar field like 
pressure, the magnetic field of a sunspot cannot be kept in equilibrium simply by 
pressure balance at the surface: the tension in the magnetic-field lines continuing 
below the surface exerts forces as well. The question of spot equilibrium thus 
involves deeper layers, down to wherever the field lines continue. This is the well- 



known "anchoring" problem of sunspots (see Parker 1979): at which depth, and 



by which agent is the sunspot flux bundle kept together? 
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3.1. The Anchoring Problem 



As an answer to this problem, Parker ( 1979 1 postulates the existence of a hori- 



zontal flow at some depth below the surface, converging on and flowing through 
the magnetic flux bundle of the spot. The drag force of this flow on the field would 
prevent the bundle from fragmenting. There are some obstacles to this idea. A 
horizontal flow is observed around spots (the moat flow) , but it is of the opposite 
sign to the proposed inflow. There is no theory for the cause of the proposed 



flow of Parker (1979). In fact, the "heat flux blocking" by the sunspot would 



cause a flow of opposite sign, and this is actually observed in the form of the 
moat flow. It is also questionable if the proposed flow would actually be sufficient 
to keep the flux bundle together, in view of the (interchange) instabilities to be 
expected in this picture ( ,Schiissler, 1984 ). What has kept Parker's proposal alive, 
however, is the helioseismic inference (Duvall et al., 1996) of a huge downflow 
flow of 2 km of the sign proposed by Parker. It is a puzzling observation, 
which if confirmed, would require new theoretical ideas. However, more recent 
helioseismic inferences derived from inversions of subsurface flows around the 



sunspot in AR 9787 (see Section 7.6.2 ) do not conflrm the existence of this 
downflow. 



A second proposal (already implicit in Cowling 1953 see also Babcock 1963 
Leighton 1969 Spruit and Roberts 1983| is that the flux bundle of a spot 



actually continues all the way to the base of the convection zone, to the layer of 
toroidal field from which the active region erupted. At this depth, a field strength 
of « 80 kG is inferred from flux tube emerging calculations {e.g. D'Silva and 



Choudhuri 1993 Caligari, Moreno-Insertis, and Schiissler 1995). The Alfven 



travel time for a change at the base to propagate to the surface (more accurately: 
the propagation time of a transverse tube wave) is then about five days. This is 
the time scale on which the spot would change if there were nothing to maintain 
conditions at the base. The idea is thus that a spot, while in equilibrium at 
the surface, would actually be a transient structure in the layers from which its 
magnetic field originates. This proposal agrees roughly with the lifetime of most 
small spots. It is too short, however, to explain the anchoring of stable long-lived 
spots. 



3.2. The Issue of Flux Emergence 



The anchoring problem of long-lived spots is thus still somewhat open. Quite 
clear, however, is the picture of flux emergence: the process by which a loop of 
flux rises from a horizontal layer of magnetic flux at the base of the convection 



zone to form an active region, as discussed briefly in Section 2.1.1 Simplified ID 
calculations in "thin tube" approximation (D'Silva and Choudhuri, 19931 Fan, 



Fisher, and McClymont, 1994 Caligari, Moreno-Insertis, and Schiissler, 1995) 
point to a field strength of about 10'' G at the base of the convection zone. 
At this strength agreement is reached with three independent key properties of 
active regions: i) the time scale of emergence (a few days), ii) the heliographic 



latitude range of emergence, and in) the tilt of active region axes {e.g. Caligari, 



Schiissler, and Moreno-Insertis 1998 and references therein). This is also the 



SOLA: SDLA1202R2_Moradi.tex; 24 August 2010; 1:25; p. 16 



Modeling the Subsurface Structure of Sunspots 



field strength at which the horizontal field is first expected to become unstable 
to the development of bends in the field lines, creating loops rising to the surface 



(Schiissler et ai, 1994). 



However, the picture presented by rising flux tube simulations is somewhat 
incomplete, since the last stages of the flux emergence process in the near pho- 
tospheric layers are typically excluded. It has been pointed out by [Schiissler and 
[Remp el, (2005) that the typical size of active regions correspond to a wavenumber 
of m = 10 to 60, while buoyant flux tubes typically prefer m = 1 or 2 instabilities. 
From this mismatch in scales one would expect an unrealistic drift of sunspots 
apart from each other much further than actually observed due to magnetic 
tension forces. Possible solutions could include a "dynamical disconnection" as 



suggested by Schiissler and Rempel (2005) or a flux emergence process that 



prefers larger wave numbers from the beginning. Observational evidence for this 



process has recently been provided by Svanda, Klvaiia, and Sobotka ( 2009 ) 



Overall, the convergence of these lines of evidence supports the basic cor- 
rectness of the flux emergence picture of active region formation that has long 



been intuitively evident from observations (Cowling, 1953). It has the important 
implication that the energy density of the magnetic held at its source is more than 
two orders of magnitude larger than energy density in convective flows; at such a 
strength the field is nearly unaffected by convective flows. The same is the case at 
the surface of a sunspot. This is incompatible with the cornerstone of traditional 
mean-fleld dynamo models based on the effect of convective turbulence acting 
on magnetic fields. The picture sketched above and the mean field convective 
dynamo model thus cannot both be correct. This fact is remarkably politely 
ignored in discussions on theories of the solar cycle. 

3.3. Anchoring Where? 

Some details of the evolution of an active region give additional clues about the 
anchoring of sunspots. After their formation, long-lived spots wander about a 



Gizon et al. 2009a 


2009b 


1, 


few days {e.g. 


Mazzucconi, 



Coveri, and Godoli[ 1990). This is similar to the inferred Alfven travel time 
from the base of the convection zone and the surface. The observed settling 
process thus agrees with the notion of anchoring in deep layers. The "settling 
time scale", the Alfven travel time from the anchor to the surface, agrees with 
the field strength at the base of the convection zone inferred from the other 
properties of the emergence process mentioned above. 

While anchoring at the base of the convection zone thus agrees with a range 
of observational indications, this does not in itself prove that other locations 
within the convection zone are excluded. Such locations are, however, somewhat 
unlikely. Apart from the boundary with the stably stratified interior it is hard 
to find a plausible location for anchoring a magnetic field in a convecting stellar 
envelope, which is itself unstable to fluid displacements. The anchor of a long- 
lived spot needs to act over a time longer than the convective turnover time of 
the envelope as a whole (which would be on the order of a month, in the classical 
mixing length view). 
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4. Sunspot Models 

4.1. Semi-Empirical Models of the Sunspot Atmosphere 

Semi-empirical models of umbral or penumbral atmospheres give the variation 
of thermodynamic variables and magnetic-field vector with optical depth based 
on empirical data and theoretical considerations of mechanical equilibrium and 
radiative transfer. Semi-empirical models can be divided into two groups based 
on the observational material and methods used for their derivation. 

One group of such models are spatially unresolved one-dimensional models 
based on multi-wavelength observations of the continuum, weak and strong 
spectral lines. The general procedure in constructing a model atmosphere usually 
involves first determining a temperature-optical depth relation, as a best fit to 
the empirical data, and then to determine the gas and electron pressures by 
integrating the equation of hydrostatic equilibrium. In general, a number of as- 
sumptions/generalizations are made when calculating this type of semi-empirical 
models: i) For the lower photospheric layers, measurements of the center-to-limb 
variation of continuum intensity at several spectral intervals and the profiles of 
the weak spectral lines (or the wings of strong lines) are used together with the 
LTE assumption. While for the high chromospheric layers strong spectral lines 
are used like Can H and K and IR lines. Ha, Nai D, etc. and the analysis is 
carried in NLTE; ii) Hydrostatic equilibrium is assumed. The question of the 
magnetic-field distribution is not addressed. 

Most models consist of a single-component, meant to represent a horizontal 
average over the umbra or penumbra, while a number of models consist of two- 
components, designed to treat bright and dark regions separately in observations 
that do not fully resolve them spatially. 

Another group of models uses a relatively new strategy based on the inversion 



of the radiative transfer equation of polarized light {e.g. Ruiz Cobo and del Toro 



Iniesta 1992). The input material for these models are spectro-polarimetric or 
spectral observations (typically high-resolution). Strict mathematical methods 
are used to find the best fit to the spectral line profiles and only few spectral lines 
are required to produce a model atmosphere. This method provides sunspots 
models (umbra, penumbra and the surrounding plage) spatially resolved to the 
extent allowed by the observations. The following assumptions are typical: 

• The inversion of radiative transfer equation is done under LTE conditions 
for photospheric lines. NLTE inversion results have been reported for the 



chromospheric ionCaii infrared triplet (Socas-Navarro, 2005). All variables 
are derived as a function of optical depth. 

The gas stratification is obtained at heights of formation of spectral lines, 
at most sal20 km below, and up to R:!300-400 km above the continuum 
level. The data about the chromospheric layers are less frequent (see e.g. 



Socas-Navarro 2007) 



The models include magnetic-field vector and line of sight velocity compo- 
nent with or without gradients. 

For practical reasons, the atmosphere is assumed to be in MHS equilibrium 
at each spatial location. 
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• The Wilson depression is not an inherent part of the models, but instead 
must be determined from additional considerations (e.g. lateral pressure 
balance) . 

Such semi-empirical models can include one, two or more magnetic components 
describing the variations of thermodynamic parameters, velocities and magnetic- 
field vector at each individual spatial location, supposing these variations are 
not completely resolved by observations. The shape of the polarization profiles 
in the penumbra suggest that they result from at least two different magnetic 



components, one of them carrying the Evershed flow (Bellot Rubio, 20031. Note 



that the magnetic-field strength and its inclination in the penumbra are cor- 



related with the brightness of the penumbral filaments ( 


Beckers and Schroter, 


1969 


Schmidt et al, 1992 


Title et al, 1992 Bellot Rubio, 2003 


Solanki, 2003 


1, 



thus the components with different field strength also typically have different 
temperatures. 

All these models are useful in constraining certain physical processes that 
determine the structure of a real sunspot atmosphere, and also in providing a 
background model for studies of element abundances, wave propagation, and 
other behavior in a sunspot. 



4.1.1. Umbral Models 



Below is a summary of some of the more often used umbral models. The first 
three belong to the first group of models described above and only provide 
variations of thermodynamical quantities: 



Avrett ( 1981 ): Provided a combined model of the umbral photosphere, chro- 



mosphere and transition zone, by combining the low photospheric model 



of Albregtsen and Maltby ( 1981 ) with the upper photospheric and chro- 



mospheric parts of Lites and Skumanich ( 1982 ) model and the transition 



region model of Nicolas et al. (1981) 



Staude et al. (1983): Derived a comprehensive umbral model covering the 



full height range from the photosphere to the corona. Their model was 
based on a large set of observations (radio, optical, EUV, X-ray). The pho- 



tospheric structure was taken largely from Stellmacher and Wiehr (1975), 



while its lower chromospheric stratification is similar to that of Teplitskaya,| 



(hot and cold) in higher layers. 



Grigor'eva, and Skochilov (1977). They also introduced two components 



Maltby et al. (1986): Came up with an improvement of the Avrett (1981) 



model in the deep layers. These consist of a set of three models (each for 



a different phase of the solar cycle). Modifications of the Maltby et al. 
(1986) models have been proposed by Caccin, Gomez, and Severino ( |1993 ), 



Severino, Gomez, and Caccin ( 1994[ ) and Ayres '{ 1996 1 , including constraints 



from more spectral lines. 



Collados et al. (1994): Presented models of the umbra of large and small 
spots from the inversions of Stokes / and V spectra from the darkest cores 
of three different sunspots. The run-with optical depth of temperature, 
magnetic-field vector, and velocity along the line of sight were obtained. 
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Their observations confirmed that there are noticeable differences in the 
temperature and field strength between the umbra of large and small spots. 

The temperature stratifications of all of the above-mentioned umbral models 
reveal an umbra which is cooler than the surrounding field-free regions in the 
photospheric layers, but slightly hotter in the chromosphere, while the transition 



to coronal temperatures takes place at a lower height in the umbra (see Staude 
1986 Solanki 1990 and Collados et al. 1994 for comparison plots). 



4-1-2. Penumbral Models and Spatially Resolved Models 

Below is a brief summary of some semi-empirical average penumbral models 
and models obtained from the inversion of the high-resolution polarized spectra 
including complete sunspots or their parts: 



Makita and Morimoto (1960), Kjeldseth Moe and Malt by (19691, Yun, 



Beebe, and Baggett (1984) and Ding and Fang (1989): Produced one 



component models, with the latter two being primarily aimed at modeling 



the average penumbral chromosphere. Kjeldseth Moe and Maltby (1974) 



produced a two-component model of similar kind that provided temperature 
as function of optical depth in dark and bright penumbral filaments. 



del Toro Iniesta, Tarbell, and Ruiz Cobo (1994): Applied an inversion 



technique to a series of high-resolution filtergrams, scanning a magnetically 
insensitive Fe i line to study the vertical temperature and velocity structure 
over a small penumbral region. They provided a mean penumbral model 
atmosphere as a function of optical depth. 



Westendorp Plaza et al- (2001): Besides deriving the velocity and tempera- 



ture profiles, they also derived the magnetic- field stratification of a sunspot 
from inversions of observations of magnetically sensitive Fe i lines at 630 nm 
with a spatial resolution of approximately 1" and one magnetic component 
together with stray light component in each observed pixel. Similar models 



Fei spectral lines at 1.56 /Ltm. Mathew et al. 



were provided by Mathew et al. ( 2003 2004 ) from the inversion of infrared 



( 2004 1 also calculated maps of 



Wilson depression and plasma parameter /3. 



Rouppe van der Voort (2002): The observed radiation temperatures in 



the Can K wing were used to derive the temperature stratification of 
fine-structure elements in the penumbra. Three atmospheric models were 
constructed to represent cool, intermediate, and hot features within the 
penumbra, with temperature differences of order 300 K between them. 



Bellot Rubio (2003) and Borrero et al. (2004 2005 20061: Inverted the po 



larization profiles from the penumbra and considered the uncombed penum- 
bral model with two magnetic components, where the component with a 
more horizontal and weaker field harbors the Evershed flow. 



Socas-Navarro ( 2005 1 : Produced a 3D sunspot model, up to chromospheric 



heights, from NLTE inversion of the Ca ii infrared triplet polarized spectra. 



Puschmann, Ruiz Cobo, and Martinez Pillet ( 2008 ) : By assuming that the 



fine structure of the penumbra is spatially resolved in high-resolution Hin- 
ode/SOT data, they invert the data with a single component model. The 
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advantage of their modeling is a very detailed and self-consistent translation 
from optical depth to geometrical height scale taking into account terms 
of the Lorentz force and using a generic algorithm. The Wilson depression 
was derived this way and the model was checked to be divergence-free and 
in equilibrium in the horizontal and vertical directions. 

Unlike quiet Sun models, sunspot models should not be used indiscriminately 
to all cases and can only be taken as representative examples. As was discussed in 
Section ["2.3[ a number of observations have shown that the thermal stratification 
of sunspots depends sensitively on its magnetic-field strength, sunspot size, and 
possibly, on the solar cycle as well. 

4.2. MHS Models 

The simplest static models of a pore or a sunspot ignore azimuthal variations 
and treat it as an axisymmetric, poloidal magnetic field (with no azimuthal 
component) confined to a homogeneous flux tube of circular cross-section. The 
equation (in cgs units) describing the MHS equilibrium of the flux tube is: 

-Vp + /7g+ ^(JxB) = 0, (1) 
c 

where p is the gas pressure, p is the gas density, g is the acceleration due to 
gravity, c is the speed of light, and B is the magnetic-field vector. The electric 
current density is given by Ampere's Law: 

J = £(VxB). (2) 

In reality, magnetic fiux concentrations are devoid of symmetry (like most 
sunspots). However, occasional long-lived spots, sufficiently separated from other 
large flux concentrations, tend to be round and to have regular, ring-like penum- 
brae. Therefore, an isolated, axially symmetric spot with a unipolar magnetic 
field is somewhat justified. Furthermore, the lifetime of sunspots is much longer 
than any dynamical timescale in the solar photosphere. Some sunspots persist for 
several rotations, essentially unchanged, whereas perturbations that propagate 
with the Alfven speed would need about an hour to cross a large sunspot at its 
photosphere. 

There are also a number of other assumptions/generalizations to keep in 
mind when considering MHS models: i) The MHS models are generally limited 
to two dimensions, the vertical and radial directions; ii) Dynamic phenomena 
(important for shaping small-scale magnetic and thermal structure) and all 
fiuctuations related to convective motions are usually ignored; in) Most MHS 
models (unrealistically) tend to treat the force balance in isolation of the energy 
balance; iv) Unless the energy equation is solved together with the force balance 
equation, either the magnetic field, temperature, or gas pressure must be spec- 
ified, as well as an equation of state; finally, v) For models in which the force 
balance and energy equations are consistently solved throughout the sunspot the 
convective transport is treated by applying the mixing length formalism and the 



radiative transport by the diffusion approximation {e.g. Chitre 1963 Deinzer 



1965 


Chitre and Shaviv 


1967 


Yun 


1970 Jahn 



19891 
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4-2.1. Force- Free and Potential Field Models 



These are the simplest static model where the magnetic field within the sunspot 
flux tube is assumed to be force-free, resulting in the atmosphere being hori- 
zontally stratified with variations only in the vertical direction. A number of 
conditions must be satisfied in order to create such a model: i) the force- free 
field must satisfy the condition: 



(V X B) X B = 0, 



(3) 



and a) the assumed axial-symmetry also requires that the field satisfy the 
current-free condition V x B = 0, and hence B is a potential field B = — V0 
(where (j) satisfies Laplace's equation, V^(/) = 0). 

Potential field models can thus be effortlessly derived by taking solutions 
of Laplace's equations {e.g. using dipole or Bessel function potentials). Simple 
models of pores were constructed Simon and Weiss ( 1970), using Bessel function 



solutions of Laplace's equation. Subsequent advances of this model were made 



by Spruit (1976) and by Simon, Weiss, and Nye (1983), who represented the 



field in a pore by a potential field such that, at the photosphere, Bz was uniform 
over a disc with a prescribed radius and zero outside it. 

One obvious advantage of this method is that direct solution of equation 
([!]) is avoided, however the solution now involves the solution of a non-linear 
boundary value and free-surface problem (i.e., determining the location of the 
current sheet, as potential fields fill the whole atmosphere unless bounded by a 
current sheet), which is nontrivial. 



.2.2. Current Sheet Models 



A more satisfactory potential field model can be derived by constructing a field 
contained within a flux tube such that the difference between the internal [pi] 
and external [pc] gas pressures is balanced by the Lorentz force in a thin current 
sheet bounding the flux tube, such that 



Pi + 



8^ 



Pe 



(4) 



This bounding current sheet is referred to as the "magnetopause" . 

Approximate solutions have been found and applied to sunspots and pores by 



Simon and Weiss (1970) and Simon, Weiss, and Nye (1983). Wegmann (1981) 



produced the first general solution to the free boundary problem and Schmidt 



and Wegmann ( 1983 1 were the first to apply this technique to sunspots, suc- 



cessfully modeling pores. More advanced current sheet models have succeeded 
in producing relatively realistic sunspot models, including the penumbra. These 
models are briefly summarized below: 



Jahn (1989): Provided an extension of the Schmidt and Wegmann (1983) 



model by introducing body/volume currents (distributed in the outer parts 
of a flux tube, below the photosphere of the penumbra), in addition to 
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a current sheet at the sunspot-quiet Sun boundary. The body currents 
contribute to the lateral force balance and affect the pressure stratification, 
so that the gas in the penumbra is hotter, thus layers of equal gas pressure 
assume higher levels than in an umbra where the field is current free. A fit 



to the magnetic and photometric profiles (taken from Beckers and Schroter 



1969 ) provided a distribution of the electric currents in their model. 
Jahn and Schmidt^ ( .1994) : Introduced two current sheets, one at the outer 
boundary of the flux tube (the magnetopause) and the other at the interface 
between the penumbra and the umbra. They did this in order to obtain a 
more realistic thermal structure of the sunspot with distinctively different 
umbral and penumbral thermal mechanisms. They also assumed that the 
umbra is thermally insulated from the penumbra, while some of the energy 
radiated from the penumbra itself is supplied by convective processes that 
transfer energy across the magnetopause. 

Pizzo ( 1990D: Extended the earlier work of Pizzo 1986 (see Section 4.2.4), 



using multi-grid techniques to calculate the magnetic structure of a sunspot 
bounded by a current sheet. 

In general, current sheet models are consistent with the concept that sunspots 



represent discrete, erupted, magnetic entities (Solanki, 2003). Observational sup- 
port for the current-sheet description of spots is taken from the sharp transition 
between the umbra and quiet photosphere in pores and from the relatively 



uniform photometric appearance of most umbrae (Gokhale and Zwaan, 1972). 
The fact that the magnetic field is so large at the white-light boundary of the 
sunspot also strongly suggests that sunspots are bounded by a current sheet 



(Solanki and Schmidt, 1993). However, as Solanki (2003) points out, the rugged 
nature of the sunspot boundary in white-light images means that the current 
sheet is not as well defined as one might picture on the basis of simple flux-tube 
models. Further evidence for a current sheet surrounding sunspots comes from 
observations suggesting that the field inside sunspots is close to potential. This 
suggests that the currents bounding the strong field must be mainly located in 



a relatively thin sheet at the magnetopause {e.g. Lites and Skumanich 1990). 
4.2.3. Self-Similar Fields 



There is considerable arbitrariness in assigning the distribution of the volume 
current and the corresponding radial structure in the sunspot atmosphere. As 



first shown by 


Schliiter and Temesvary| ( 


1958|) (and later extended by others. 


e.g. 


Chitrel|1963| 


Jakimiec 1965 


Jakimiec and Zabza| |1966| |Chitre and Shaviv| 


1967 


), the problem can be greatly simplified by assuming a self-similar profile 



for the magnetic field. 

In cylindrical coordinates, the (untwisted) Bz and Br components of the 
magnetic field take the form: 



BAr,z) = .fiC)Bo{z) 



(5) 



Brir, z) 



2^^^' dz ■ 



(6) 
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where r and z refer to the radial and vertical coordinates respectively, Bq{z) 
is the field strength at the flux-tube axis, and C = f ■\/-Bo(^)- The shape of the 
function /(C) may be freely chosen (usually a Gaussian). For a non-constant 
f{C), inserting Equations ^ and ^ into ([T]) reduces the equation of MHS 
equilibrium to the following equations: 

dp fdBr dB, 



dr An \ dz dr 



(7) 



dp Br fdBr dBA 

Integrating Equation Q over r from to infinity for constant z leads to the 
following expression for Ap(z) = Pe{z) — Pi{z), the difference in gas pressure 
between the external [p{oo,z)] and internal [p(0, z)] regions: 

Svr \27r az^ J 

where 4> denotes the total magnetic flux and y — \J Bo{z). Thus this method 
has the advantage in that it simplifies the mathematical treatment of MHS equi- 
librium by reducing the partial differential equation to a second-order ordinary 
differential equation for the field strength at axis of the spot by specifying the 
cross-sectional shape of the magnetic field distribution within an axisymmetric 
flux tube. 

A general assumption made is that the distribution of magnetic flux on hor- 
izontal planes is geometrically similar at each depth. Furthermore, in contrast 
to current sheet models, self-similarity allows for a continuous variation of field 
strength and gas pressure across the spot. The field falls off smoothly from the 
central axis value to zero at large radial distances (hence, no clear definable 
"inside" or "outside" of the spot in this description, therefore essentially ignor- 
ing the fact that sunspots have sharp edges). Thus, there is the computational 
convenience as treatment of the discontinuity associated with a current sheet is 
avoided. However, the similarity law enforces a somewhat arbitrary distribution 
of electric current, resulting in the appearance of a bright ring in the emergent 
intensity. Those currents determine (to some extent) the horizontal temperature 
variations, which is in general need not comply with the observed photometric 



profile (Jahn, 1992) 



A further disadvantage of a purely self-similar sunspot model is that nega- 
tive pressures and densities are often obtained in the photospheric and upper 
atmospheric layers of the flux tube, due to the fact that the hydrodynamic 
pressure and density are decreasing exponentially with height, while the mag- 
netic field does not quite decrease at the same rate. There have been some 



proposed workarounds to this problem (see e.g. Hanasoge 2008), however such 
methods are not ideal as they tend to substantially alter the governing differential 
equations because they require the inclusion of terms in the ideal MHD equations 
which are not physical. 
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Furthermore, since the Lorentz force also drops with height, the magnetic 
field essentially becomes an unbounded, force-free field within a few pressure 
scale heights above the photosphere, resulting in a field configuration not too 
different from a potential field. This fact implies that self-similar models are 
essentially not force-free in the regions where they should be. Since the shape of 
the magnetic field at the surface is sensitive to this, one has here a model that 
breaks down in an essential way in just those regions where diagnostics are best. 

Nonetheless, due to their simplicity, similarity expansions have been utilized 
to generate the field configuration for a number of studies, some of which are 
summarized below: 



Deinzer (1965): Generalized the Schliiter and Temesvary (1958) model. 



where the similarity law for the magnetic field is coupled with the thermo- 
dynamic structure along the axis of a spot, as described by mixing length 
theory. 

Yun ( 1970j ): Improved the Deinzer (1965) model by the introduction of an 



"effective surface monopole" , which controlled the inclination of the partic- 
ular field lines identified with the outer edge a spot at the surface. Hence, 
the upper boundary condition takes into account the fact that the gas 
pressure difference at the photospheric level is not negligible (as assumed 



by Deinzer 1965 ) . Lower boundary conditions were also modified, as the 
effects of partial ionization on the relation between the internal and external 
pressures and temperatures were included. 



Yun (1971) and Osherovich and Flaa (1983): Demonstrated that the intro- 



duction of a moderately twisted field (« 17° near the surface, compatible 
with observations) contributes little to the force balance in spots and only 
slightly changes the main characteristics of the model (as already mentioned 
in Section 



2.2.2). 



Landman and Finn (1979): Imposed an Evershed-type radial velocity dis- 
tribution in the upper region of the spot atmosphere, in order to get a 
satisfactory continuum intensity profile across the sunspot. However, rel- 
atively large values of the Evershed flow (i.e., close to 10 km s~^) were 
required to obtain a satisfactory temperature profile. 



Low (19801: Prescribed a method for generating exact solutions of MHS 
equilibrium describing a cylindrically symmetric magnetic flux tube ori- 
ented vertically in a stratified medium. Given the geometric shape of the 
field lines, compact formulas were presented for the direct calculation of all 
the possible distributions of pressure, density, temperature and magnetic- 
field strength compatible with these field lines under the condition of static 
equilibrium. A particular solution was obtained by this method for a mediumj 



Temesvary 


(1958 


Osherovich 


(1982 



Osherovich ( 1982 1 : Extended self similar models to include field lines in the 
outer part of the sunspot that return to the solar surface just outside the 
visible sunspot (return flux). The emerging flux constitutes the penumbra. 
The predicted continuum intensity of return flux models was not much 
closer to the observations than the standard self-similar models. 
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Fla, Skumanich, and Osherovich ( 1982 ) : Applied the return-flux model to 



a spot with the observational data of Lites and Skumanich ( 1982 ) for 



pressure, maximum field strength and size. The force balance equation 
was solved to obtain self-consistent magnetic field, pressure, and tempera- 
ture distributions. The resulting distributions appeared to yield improved 
representations of unibral-penumbra and penumbra-quiet Sun boundaries 



compared to regular {e.g. Schliiter and Temesvary 1958) self-similar mod 



els. However, it appears that one needs to introduce an Evershed flow to 
eliminate the apparent umbral bright ring in the continuum. Similar work 
on return-flux sunspot models was undertaken by Osherovich and Lawrenc"e] 
(|1983|),|Osherovich and Garcial (|1989|) andlLiu and Songl (|1996|). 



Solov'ev (1997): Extended the self-similar sunspot models by introducing a 



current sheet at the sunspot boundary. 



Moon, Yun, and Park (19981: Included a description of the energy balance 



and the observed horizontal variation of the Wilson depression when de- 
termining the shape function from the observed radial dependence of the 
magnetic field 



Cameron, Gizon, and Duvall (2008 



and Gaily (20091, and Shelyag et al 



Hanasoge (2008), Moradi, Hanasoge, 



( 2009 1 : All employed simple self-similar 



toy sunspot models in conducting numerical simulations of helioseismic 
wave propagation through magnetized plasmas. 

As we shall see in the proceeding sections however, constructing more physically 
realistic sunspot models is practical, especially with the MHD codes (Sectionjoj) 
that are currently used for the wave-propagation problem. 

4.2.4- Solution of Full MHS Force Balance 



These methods involve solving for the magnetic field on the basis of full MHD 
equilibrium, with and without a current sheet. Pressure is specified throughout 
the numerical domain (usually as a function of depth and taken from semi- 
empirical models) , partly depending on the distribution of field lines (hydrostatic 
equilibrium acts along each field line). 



Pizzo ( 1986) utilizes the description of Low (1975), who proposed transform- 



ing the thermodynamic parameters (the pressure and temperature of the gas) 
into functions of the magnetic vector potential and depth. In this form, a func- 
tional form is prescribed for the gas pressure, but not for the magnetic field as 
in the similarity models, and the equilibrium is solved as a classical non-linear 
boundary value problem. 

By transforming the pressure and density into functions of the field line con- 
stant u (used in both |Low[ |1980 and Pizzo 1986 it essentially determines the 
shape of the field lines) and height z, and requiring hydrostatic equilibrium along 
the field lines. Equation [T] can be reduced to a single scalar equations describing 
the magnetostatic equilibrium of an axisymmetric, poloidal field: 



o u 

Qj.2 



1 du 
r dr 



o u 



: dP{u,z) 
du ' 



(10) 
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where P is a function related to the gas pressure distribution. Low ( 1975 1 pro- 



vides an approximation for the distribution of gas pressure along the magnetic- 
field lines in a vertical, axisymmetric flux tube in magnetostatic equilibrium, 



P{u, z) = Pq{u) exp 



dz' 
h{i^,z') 



(11) 



where Po(u) is the gas pressure along the lower boundary, h{u^ z) is the isother- 
mal scale height {h = TZT/jig, where TZ denotes the ideal gas constant, T the 
temperature, the mean molecular weight, and g is the acceleration due to 
gravity) for a plasma obeying the ideal gas law (p = pTZT/fi). The u = con- 
stant curves describe the field lines of the system. The vertical and radial field 
components may then be expressed in terms of u: 



and 



Br 



1 du 
r dr 

1 du 
r dz 



(12) 



(13) 



The range of validity is determined by the representative pressure distributions 
along the axis and in the field-free atmosphere. 

The gas pressure difference between the quiet photosphere and the axis of 



the spot is needed for the computation. Pizzo ( 1986 1 takes these values from the 



semi-empirical models of the umbral photosphere derived by Avrett ( 1981 ). Pizzo 



( 1986| ) then develops a method for the iterative numerical solution of Equation 
(10), essentially a second order, non-linear, elliptic partial differential equation. 



which can be easily solved using standard numerical techniques in the case of 
fixed boundary conditions. 

An advantage of this model is that the Wilson depression and net internal- 
external pressure difference can be adjusted by vertical translation of the abso- 
lute height scales of the two reference atmospheres. However, the configuration 



considered by Pizzo ( 1986 1 has its base placed 120 km below the visible surface 



of the umbra, which corresponds to z = in the Avrett ( 1981 ) model, hence 



his models do not address the question of the spot structure in deeper layers. 
The model also assumes a Gaussian profile of the magnetic field across the 
base (i.e. self-similar), thus ignoring the existence of a discontinuous transition 
from the magnetized to field free plasma (i.e., a current sheet). However, Pizzo 
( 1990 ) later extends his method by incorporating the free-surface problem in the 



solution of the equation of magnetostatic equilibrium for a flux tube surrounded 
by an infinitely thin current sheet, utilizing a body- fitted mesh generation and 
multi-grid relaxation techniques for solving Equation (10 1. 



Steiner, Pneuman, and Stenflo ( 1986 ) also developed a method for the itera- 



tive numerical solution of Equation ( 10 1, including a boundary current sheet and 



also field twist in their treatment, while Gaily ( 1991 ) adopted a full multi-grid 
method to tackle the free boundary problem by formulating it in terms of inverse 
or fiux coordinates, in which the magnetic-field lines become coordinate lines. 
This results in the energy equation reducing to ordinary differential equations 
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along field lines when the radiation is optically thin. Also, if steady plasma flow is 
allowed, Alfven's theorem guarantees that there can be no cross-field component 
of velocity, i.e., that the fluid flows along field lines. 



Khomenko and Collados ( )2006) used the Pizzo (1986) sunspot model in nu- 



merical simulations of magnetoacoustic wave propagation, and more recently 



Khomenko and Collados (2008) produced a new set of models consisting of 



concatenation of self-similar models in the deep layers, where the gas pressure 
dominates over magnetic pressure, with models in which the pressure distribu- 
tion is prescribed on the axis. In the deep photospheric layers, a self-similar 



solution for the magnetic field is calculated following the method of Low ( 1980 ), 



while the pressure and density distributions with height and radius are found 
from analytical expressions. A potential solution is then generated above some 
arbitrary height using the method of Pizzo ( 1986 ) - the bottom boundary of 



this model coincides with the top boundary of the deep photospheric/self-similar 
model (approximately z = —1 Mm, can be adjusted however to limit the upper 
boundary of the self-similar model). This solution is then used as an initial guess 
in the integration of the complete force balance equation along the magnetic-field 
lines, i.e., Equation (10). 



The analytical description of the pressure distribution along magnetic-field 



lines is taken from Pizzo (1986) and Low (19751. The Model S ( Christensen- 



Dalsgaard et al., 19961 pressure distribution for the field- free atmosphere is 



smoothly joined to the VAL-C (Vernazza, Avrett, and Loeser, 1981) model of 
the solar chromosphere. At the axis, the Avrett] ( 1981 ) umbral core model is 
used in the upper layers, while the linear inversion model of |Kosovichev, Duvall,| 
and Scherrer ( 2000| ) is used for deeper layers (down to a depth of 1 Mm, thus 
ignoring the "hot" layer) , which takes the Wilson depression to be at 450 km. A 
smooth transition between the models is then calculated for the gas pressure and 



scale height distributions, in the same manner as Pizzoi ( 1986 ), and by changing 



the parameters of the solution, a set of models with the desired properties 



can be produced. However, as with the Pizzo ( 1986 ) and all other pressure- 



distributed sunspot models, the vertical extent of these models is severely limited 
by available semi-empirical data on the sunspot axis. 

4.2.5. Extrapolated Field Models 



Martens et al. ( 1996 ) present a force-free constant-a model for the magnetic field 



in and above a fluted sunspot. They demonstrate that magnetograms for round 
sunspots can essentially be matched by a series of Bessel functions. Their model 
parameters are chosen to reproduce the high resolution observations (magne- 
tograms) of Title et al. ( 1993 ) at the 1-m Swedish Solar Observatory at La Palma, 
and an analytical expression is obtained for the 3D magnetic field emanating 
from the sunspot 's umbra and penumbra. The model accurately reproduces the 
azimuthal variation in inclination angle, as well as the mean constancy of the 
magnetic- field strength, and the appearance of a highly corrugated neutral line 
on the limb side of off-center sunspots. 

Another model based on matching mnagnetograms with analytical expres- 



sions is the axisymmetric sunspot model of Moradi and Cally (2008), which 
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consists of a non-potential, untwisted, 3D MHS sunspot model constrained to 
fit observed surface magnetic-field profiles. The preferred surface field config- 
uration of the sunspot model was derived from constrained polynomial fits to 
the observed scatter plots of the radial [Br) and vertical {Bz) surface magnetic- 
field profiles of AR 9026 on 5 June 2000 - a fairly symmetrical sunspot near 
disk-center, ideal for helioseismic analysis - obtained from IVM (Imaging Vec- 
tor Magnetograph) vector magnetograms. The surface field is therefore quite 



realistic, which is important because there is evidence {e.g. Schunker and Cally 



2006) that magnetic effects in helioseismology are dominated by the top few 



hundred kilometres. The fits of Br and B^ are then used to derive an analytical 
form for the fiux function. Instead of a current sheet along the boundary, the 
authors formulate an analytical form for the outermost field line to allow field 
strength to smoothly drop to zero. However, in a similar vain to other sunspot 
models that derive their thermodynamic properties solely from a prescribed 



magnetic-field configuration {e.g. self-similar models. Section 4.2.3 1, the axial 
pressure and density of these models posses a high degree of sensitivity near the 
upper boundary, with a tendency to become negative in the photospheric layers 
and beyond. 



4.3. Non-MHS Models 



.3.1. Semi-Emprical Sunspot Models 



Very close to the solar surface, the properties of sunspots can be inferred from 
spectro-polarimetric measurements. These layers are the most important to 
model since it is exactly in this region that the effects on the wave cannot be 
treated assuming only weak perturbations. Cameron et al. (2010| have empha- 



sized the desirability of incorporating this information into helioseismic wave- 
propagation simulations by constructing simplified, axisymmetric sunspot mod- 
els which are designed to capture the effects of the first few hundred km. 

Since full spectro-polarimetric inversions are available for very few spots {e.g. 



Mathew et al. 20041, Cameron et al. (2010) utilize a combination of existing 



semi-empirical models: the umbral model of Maltby et al. { 1986 ) and the penum- 



bral model of Ding and Fang ( 1989 1 , which cover the range of heights from about 



500 km below the quiet-Sun r = 1 level to the lower chromosphere. The OPAL 
equation of state tables are used to infer sound-speed, temperature, pressure 
and density profiles. Since the aim is to only model the near-surface region. 



the atmospheric models are smoothly matched to the Model S ( Christensen- 



Dalsgaard et al., 1996) quiet-Sun atmosphere below tsqoo = and 400 km above. 



Cameron et al. (2010) combine these ID models to form a 3D, axisymmetric 



model of the sunspot in AR 9787, with the umbral and penumbral radii cho- 
sen to match those of the observed sunspot (10 and 20 Mm respectively). The 
model is then stabilized with respect to convection by calculating the relative 
perturbations with respect to the quiet-Sun atmosphere and imposing these on 



the stabilized model (see Section 6.2 for details). As good measurements of the 



magnetic field are unavailable for the sunspot in AR 9787, the model assumes 
that Bz is described by a Gaussian horizontal profile and the vertical profile is 
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Figure 3. Thermodynamic, sound and fast mode speed profiles of the sunspot model of AR 
9787 produced by Cameron et al. (2010). Left panel: Vertical temperature profile of reference 
stratification (solid) and near the center of the umbra (dashed). Right panel: Vertical sound 
speed profile of reference stratification (solid) and near the center of the umbra (dashed). The 
dotted line shows the fast mode speed in the center of the umbra. 



tuned so the field inclination at the umbra-penumbra border is approximately 
45°. Figure [3] shows the vertical temperature, sound and fast mode speed profiles 
in the center of the umbra of the sunspot model. 

4-3.2. Dynamical Models 



Kitchatinov and Riidiger ( 2007 ) propose a theoretical sunspot model, with com- 



plete dynamics of both magnetic field and flow, which can reproduce the observed 



bright rings around sunspots {e.g. Rast et al. 1999 2001). A simphfied model 



was initially used to probe the possibility of reproducing bright rings by heat 
transport alone. In this model the mean flow velocity is put to zero, the spot- 
like structure of the magnetic field is prescribed and stationary, and the diffusion 
equation for the entropy is solved. The authors find the brightness of the region 
occupied by magnetic field decreases due to magnetic quenching of the thermal 
diffusivity, however the resulting surface brightness profile did not show a bright 
ring. 

A more consistent MHD model of a sunspot was then considered using the 
complete momentum equation together with the induction equation. Amplitudes 
of the initial uniform field of several hundred Gauss were considered. The surface 
flow generated is convergent near the spot but divergent at larger distances. The 
amplitude of the flow modeled was approximately 800 m s^^. Both brightness 
and fleld strength (central value of about 2700 G) were modeled to be almost 
uniform in the central parts of the sunspot and changing rapidly with radius be- 
yond. The radial heat-flux profiles from their simulations show consistent bright 
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rings around the spots, appearing to be somewhat brighter than observed rings. 
To probe the contribution of the flow to the bright rings, the field is switched 
off. The bright rings do not disappear, which leads the authors to conclude that 
both the flow and the reduced diffusivity quenching contribute to the resulting 
bright rings. 

Another set of dynamical MHD sunspot models, consisting of idealized ax- 
isymmeric flux tubes (where the magnetic field is matched to a potential field at 
the upper boundary) in a compressible convecting atmosphere, were presented 



by Hurlburt and Rucklidge (2000). In their models, they find the magnetic flux 
to be confined by an inward "collar" flow at the surface. Further outside, the flow 
direction appears to reverse and a moat cell appears. The authors suggest that 
the collar cell holding sunspots together is hidden beneath their penumbrae, 
so that only the outflow in the moat cell is visible at the surface. Although 
this particular flow pattern (i.e., inflows and down- flows around sunspots, first 
proposed by |Parke"r} 1979), can in theory help to explain the question of the sta- 
bility of long-lived sunspots, there are a number of theoretical concerns with this 
conjecture (already discussed in Section|3|. Furthermore, the existence of such 
a flow structure around sunspots is absent from both observations and recent 



simulations of sunspot structure (Heinemann et al. 20071 Rempel, Schiissler, 



and Knolker 2009a Rempel et al. 



helioseismic inferences (see e.g. Section 7.6.2 ), as well as realistic radiative MHD 



2009b see below for details). 



4.4. Numerical Simulations of Radiative Magnetoconvection 



Significant progress in our ability to simulate sunspots using realistic MHD 
simulations {i.e. MHD simulations that include the solar equation of state and 
multi-dimensional radiative transfer) was only possible over the past couple 
years. This is primarily due to the fact that pursuing radiative MHD simulations 
on the scale of sunspots with sufficient resolution for capturing the essential 
scales of magneto-convective energy transport requires fairly large computational 
domains and accordingly computing power. Additional to the computational 
domain size also the physical parameters encountered in and above the umbral 
region of a sunspot pose significant numerical challenges. The combination of 
several kG magnetic field with the rather small density scale height leads to a 
steep increase of the Alfven velocity above the sunspot umbra reaching values 
in excess of a few 1000 km s~^. Such high velocities lead to severe time step 
constraints for explicit codes that make such a simulation almost impractical 
unless this constraint is relaxed by artificially limiting the Lorentz force in low 
beta regions. 

3D models with a slender rectangular slab geometry have been used to study 



the inner penumbra (Heinemann et al., 2007 Rempel, Schiissler, and Knolker, 



2009a). The results from the MHD simulations show the formation of filamentary 



structures resembling those in the inner penumbra of a real sunspot, including 
bright filaments containing central dark cores. More recent 3D calculations of 



round sunspots, and even pair of opposite polarity spots (e.g. Rempel et al 



2009b ) , also show extended outer penumbrae with realistic Evershed flows having 
mean flow velocities of up to 6 km s^^. 
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Figure 4. Simulated sunspot in a 50 X 50 X 8 Mm domain. The top panel displays the 
bolometric intensity, the bottom panel field strength on a cut through the center of the spot. 



Most of these simulations are focused around the sunspot fine structure (um- 
bral dots and penumbral filaments) and cover only a temporal evolution of a 
few hours. Simulations addressing the long term evolution of sunspots (on a 
time scale of days) are currently only feasible if the resolution is decreased, 
therefore not resolving details of the fine structure as well as penumbral regions. 
Figure |4] displays the bolometric intensity and subsurface magnetic-field strength 
for a sunspot evolved over a total of 15 hours. Figure [5] compares the profiles of 
temperature, sound speed and fast mode speed in the center of the spot to those 
of the reference stratification near the edges of the domain. The temperature 
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profile (left panel) shows a Wilson depression of about 600 km which forms 
rather quickly during the first few hours of the simulation, deeper down a slowly 
progressing cooling front leads to a moderate adjustment of temperature in layers 
1-2 Mm beneath the t — 1 level in the umbra corresponding to changes of the 
speed of sound (right panel) on the order of a few 100 m s^^ to about 1 km s""'^. 
The fast mode speed shows a steep increase above the r = 1 level in the umbra, 
the maximum speed is here artificially limited to about 60 km to relax the 
stringent Alfvenic time step constraint. 





Figure 5. Left panel: Vertical temperature profile of the reference stratification (solid) and 
near the center of the umbra (dashed) of the sunspot model in Figure[4] The large temperature 
perturbation corresponds to a Wilson depression of about 600 km. Right panel: Vertical sound 
speed profile of reference stratification (solid) and near the center of the umbra (dashed). The 
dotted line shows the fast mode speed in the center of the umbra. 



While these simulations provide a realistic description of the thermal and 
magnetic structure in the upper most few Mm of a sunspot, they cannot address 
fundamental questions regarding the subsurface structure and origin of sunspots. 
Currently the domains are with 8 Mm still rather shallow (deeper domains up 
to 16 Mm are in preparation) and the initial field structure is based on a mono- 
lithic model. Furthermore the magnetic field is fixed at the bottom boundary. 
Nevertheless these models provide very useful artificial oscillation data that can 
be used to test helioseismic inversion methods, since the latter requires primarily 
a consistent rather than a fully realistic model. 



5. Diagnostics Potential of Helioseismology 

In the previous section we reviewed models of sunspots; each of these provides 
a prescription for computing the subsurface structure {i.e., the 3D forms of the 
sound speed, density, and magnetic field) of a model sunspot over some range of 
depths. These prescriptions are based on a wide variety of physical assumptions. 
In general, it is not known which of these assumptions are good approximations 
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(in the sense of accurately describing solar conditions). Hclioseismology promises 
the ability to infer the subsurface structure of sunspots and active regions, and 
thus to test the physical approximations that are used in building models of 
these regions. 



Hclioseismology has shown (see Gizon and Birch 2005 for a review) that 



wave propagation is different in sunspots than in the quiet Sun. For example, 
single-bounce wave travel times are altered by (order of magnitude) about thirty 
seconds (the details depend strongly on the data analysis filters, the first bounce 
distance, and the wave frequency, see e.g. Braun and Birch 2008). In Section]?] 
we provide a follow up to the detailed helioseismic study undertaken by [Gizon] 



et al. { 2009a 2009b ) of the active region AR 9787. These measurements clearly 
demonstrate that hclioseismology can measure the changes in wave propagation 
that are associated with sunspots and active regions. 

There have been many conclusions drawn about sunspot and active region 



structure from local hclioseismology {e.g. among others. Fan, Braun, and Chou 



|1995[ IKosovichev} |1996[ [Jensen, Jacobsen, and Christensen-Dalsgaard[ |1998[ 
|Chou, Sun, and Chang"! [2000{ [Kosovichev, Duvall, and Scherrer[ [2000] [Zhao] 



and Kosovichev[ 2003 Couvidat et al. 2004 Basu, Antia, and Bogart 2004 



Crouch et al. 2005 Bogart et al. 2008). Despite the long history of work on 



this topic, there is no general agreement on the subsurface structure of active 
regions. Here we will investigate some of the possible causes of the disagreement 



shown by Gizon et al. (2009a 2009b): the frequency dependence of travel times 
(Section 7.2 ), the effects of phase-speed filters on measurements (Section [7.3[ ), 
and the effects of inversion parameter choices and inclusion of surface terms on 



ring-diagram structure inversions (Section 7.4 ). 

Another approach to using hclioseismology to constrain sunspot models is the 
direct simulation of wave propagation through model sunspots and comparison 



the resulting wave field with the helioseismic measurements (e.g. 


Cameron, Gi- 


zon, and DaiffallahH2007| [Cameron, Gizon, and Duvall|[2008| [Moradi, Hanasoge, 


and Cally 


2009 


). This approach bypasses the need for linear forward models {i. e.. 



models for the linear sensitivity of travel-time shifts or ring-diagram frequency 
shifts to changes in subsurface structure) . The approximations made in the pro- 
cess of formulating linear forward models may be one source of the disagreement 



shown by Gizon et al. (2009a 2009b). A limitation of the direct simulation 



approach is that it is extremely computationally intensive. The examples shown 



in Section 7.8 demonstrate that numerical simulations are capable of reproducing 



many features of the helioseismic measurements (see also Cameron, Gizon, and 



Duvall, 2008). In addition, simulations promise to provide a powerful tool for 
determining the sensitivity of helioseismic measurements to small changes in the 
subsurface structure of sunspot models. This type of study is required to be 
able to make meaningful estimates of the signal-to-noise ratios required for the 
helioseismic measurements to constrain any particular aspect of sunspot models. 
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6. Numerical Forward Modeling of Waves Through Model Sunspots 

6.1. Numerical Methods 

Our best chance at constraining the interior structure of sunspots comes with 
constructing accurate numerical forward models. There are two major reasons 
for claiming this: i) the departure from the quiet Sun wave propagation physics 
being so dramatic in a sunspot, it is quite likely that the single scattering Born 
approximation {e.g. Birch, Kosovichev, and Duvall 2004 1 entirely fails {e.g. 



Gizon, Hanasoge, and Birch 2006). It is therefore imperative to solve the wave 
propagation equations in a fully 3D, highly magnetized environment. However, 
the basis of semi-analytical sunspot formalisms is significantly undermined by 
the sheer mathematical complexity of dealing with the MHD equations (not for 
the want of trying; e.g. Bogdan 1999), and hence the necessity of using purely 
numerical techniques, ii) Computational solvers are flexible and can be rather 
easily extended to include more effects, such as flows, changes on wave sources, 
etc. 

The difficulties in creating an MHD solver that produces sufficiently accu- 
rate results are not to be underestimated. A number of groups are working on 



this problem ( 


Khomenko and Collados, 2006 


Cameron, Gizon, and Daiffallah, 


2007 


Parchevsky and Kosovichev, 2007 


Hanasoge, 2008 


Shelyag et al, 2009 





The numerical techniques, and validation and verification procedures vary sig- 
nificantly from one code to another. A further complexity is introduced by an 
Alfven speed that increases exponentially with height, approaching values of 
several hundred km near the computational upper boundary. Simulating 
wave propagation in a plasma so highly magnetically dominated is remarkably 



expensive, especially in 3D. A number of simulations {e.g. 


Cameron, Gizon, and 


Daiffallah 


2007 


Cameron, Gizon, and Duvall 


2008 


Hanasoge, 2008; Moradi, 


Hanasoge, and Cally 


2009 Rempel, Schiissler, and Knolkerj |2009a,) apply an 



Alfven wave speed limiter to moderate the action of the Lorentz forces, thereby 
indirectly controlling the plasma-/?. Having described this context, three rather 
serious issues remain as yet outstanding: i) What are the spatio-temporal reso- 
lution requirements for simulating waves in an MHD environment? ii) How do 
we treat the upper magnetic boundary? and Hi) What harm does the Alfven 
speed limiting term do to the near-and far scattered wave field? 

As yet, only parts of these questions have been answered. This may be at- 
tributed to a systemic shortsighted forward modeling approach in our intensely 
observationally driven field: simulate and obtain results as quickly as possible. 
Evidently, we must attempt to overturn this trend by z) ensuring that the nu- 
merical methods are high-order and highly precise, ii) validating and verifying 
that the equations are being solved accurately (i.e., testing against a number 
of known solutions), Hi) determining and matching the spatio-temporal resolu- 
tion requirements, iv) using stable, physically-motivated well-tested boundary 
conditions, and finally v) testing the approximations used in the calculations do 
not directly affect the quality of the solution. Of course, attendant questions of 
computational efficiency must also be addressed, because the forward modeling 
approach requires exhaustive testing of sunspot models, requiring a large number 
of calculations. 
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6.2. Background Models Stabilized Against Convective Instability 

Numerical simulations of seismic wave propagation through random media must 
begin with an initial background model of the Sun. In helioseismology, it is 



common practice to begin with Model S {e.g. 


Cameron, Gizon, and Daiffallah 


2007 


Parchevsky and Kosovichev 


2007 


Khomenko et al. 


20091. However, most 



background models include convection and numerical forward models which 
simulate linear wave propagation are sensitive to this. Therefore, it is essential 
to remove such convective instabilities from the background in order to be able 
to successfully simulate linear wave propagation on the time scales (« 8 hours) 
required for computational helioseismology, as well as to ensure that convective 
modes do not swamp the signatures one is interested in analyzing. 

There are a number of ways to stabilize background models against convection 



{e.g. 


Hanasoge et al. 


2006 Cameron, Gizon, and Daiffallah 


and Kosovichev 


2007 


Shelyag, Fedun, and Erdelyi 


2008), 



2007 Parchevsky 



2008), all requiring zero 



buoyancy so that a vertically displaced packet of gas in adiabatic equilibrium 
will not continue to rise, i.e., ensuring that the Brunt-Vaisala frequency always 
remains positive, specifically N"^ /g > where g is the gravitational acceleration. 
There are a number of ways to go about this, for example the method employed 



by Parchevsky and Kosovichev (20071 ensures that, when this condition is not 
met, the value is replaced with zero, or very small values {e.g. 3 x 10^^ Mm^^). 



Shelyag, Fedun, and Erdelyi (2008) and Shelyag et al. (2009) have a slightly 
different approach, whereby they adjust the pressure and density of Model S 
using the equation of state for an ideal gas to retain a constant Fi — 5/3, with 
the additional constraint that the modified sound speed profile does not differ 



substantially from the original. Both Parchevsky and Kosovichev (2007) and 



Shelyag et al. (20091 have shown relatively solar-like power spectrums. 



Another method has been employed by Cameron, Gizon, and Daiffallah ( 2007 ) 



and Cameron, Gizon, and Duvall (2008), and undergone detailed development 



and testing in Schunker, Cameron, and Gizon ( 2010) , whereby the Model S back- 
ground is stabilized against convective instability by ensuring dzP = max{c^dzPo—^ 
e,dzPo), where specifies the partial derivative with respect to height, the 
subscript "0" indicates unperturbed Model S values, and e = 10^^ g s^^ near the 
surface and zero everywhere else. These changes to the background model alter 
the eigenmodes of the problem, and this must be taken into account. |Schunker^ 
Cameron, and Gizon (20101 have attempted to do this by examining the eigen- 



modes and ensuring that they are solar-like and have successfully demonstrated 
quite solar-like power spectrums. Other ways to remove the convective instability 
have not been explored however, and it is not clear which form of stabilization 
least affects the eigenmodes. 

6.3. Numerical Codes 

In this section we provide a brief summary of the four numerical simulation 
codes that were discussed and used during the two HELAS workshops on sunspot 
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seismology held in 2007 and 2009. These codes compute the propagation of solar 
waves through prescribed background model^ 



6.3.1. The I AC MHD Code 



The lAC MHD code, described in the works by Khomenko and Collados ( 2006 ) 



and Khomenko, Collados, and Felipe (2008), solves the non-linear MHD equa- 



tions for perturbations, written in the conservative form, using a fourth-order 
central difference scheme and advanced in time by a fourth-order Runge-Kutta 



method. In a similar manner to Stein and Nordlund ( 1998 ) and Gaunt and Korpi 



(2001), in order to damp high-frequency numerical noise on sub-grid scales, the 



physical diffusive terms in the equations of momentum and energy are replaced 
by artificial equivalents. In the induction equation, the magnetic diffusion term 
is retained, with ry being replaced by an artificial value. Depending on the sim- 
ulation. Perfectly Matched Layers (PML; e.g. Berenger 1994) are also placed 



at the boundaries that absorbs the incoming waves and prevents their spurious 
reflection and return back to the physical domain. For the best results, 10 to 15 
grid points are allocated to the PML layer. This code has been used to study 



the wave propagation and refraction in a small sunspot (Khomenko and Colla- 



dos, 2006); non-linear wave propagation, shock formation, mode conversion and 



energy transport in small-scale flux tubes with internal structure (Khomenko, 



Collados, and Felipe, 2008 ) 



Recently, the code was also used to model the propagation of helioseismic 



waves below the sub-photospheric structure of sunspots (Khomenko et al, 2009 



Khomenko and Collados, 2009). The most important results obtained for helio- 



seismic wave propagation below sunspots are the following: i) the fast magneto- 
acoustic mode represents an analog of quiet Sun p-modes modified by the pres- 
ence of magnetic field, ii) helioseismic waves below sunspots are sped up by the 
magnetic field by 20-40 seconds compared to the quiet Sun, Hi) the magnetic 
field produces a strong frequency dependence of the travel times, iv) the eikonal 
solution gives a qualitatively good approximation for the numerical solution, 
and finally, (v) the high-frequency fast mode waves are refracted in the magnet- 
ically dominated layers and inject additional energy, possibly causing the power 
increase observed in acoustic halos surrounding active regions. 



6.3.2. The SUM Code 

The Semispectral Linear MHD (SLiM) code solves the ideal linearized MHD 
equations using a spectral expansion in the horizontal directions and a two-step 
Lax-Wendroff treatment in the vertical. The code includes two absorbing layers 
at the top and the bottom of the box. In the top layer the waves are heavily 
damped and the effect of the Lorentz force is systematically reduced. Likewise 
the bottom layer damps the waves that propagate downward. The code has 



^We note that other numerical simulation codes also exist (e.g. 
|2007| [Hartlep, Miesch, and Mansour[|2008[ |. 



Parchevsky and Kosovichev 
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been tested against analytic solutions, which are described in detail in Cameron, 



Gizon, and Daiffallah| (|2007[) . 

In |Cameron, Gizon, and I)uvall| ( [20081 ) and |Gizo n et a/.'f2009a| |2009bD , SLiM 
was used to study wave propagation through a simplified monolithic model 



sunspot embedded in a stabilized quiet-Sun model atmosphere (see Section 6.2 1. 
The corresponding wave field computed with SLiM was then compared with 
MDI observations of /- and p-mode scattering by the sunspot in AR 9787. The 
comparisons were quite encouraging as the numerical simulations from SLiM 
were able to reproduce wave absorption and scattering phase shifts. The code has 
also recently been used on the sunspot model described in Section ["4.3.1[ as well 
as in simulations where the magnetic field of the sunspot is essentially "switched 
off", while maintaining the sound speed, pressure and density perturbations of 
the sunspot model. This type of simulation is made possible by the fact that the 
background does not need to be in pressure balance. This type of experiment 
allows one to disentangle the contributions of the Wilson depression and sound 
speed/thermal changes from mode conversion and other magnetic-field effects. 

6.3.3. The SPARC Code 

The Seismic Propagation through Active Regions and Convection (SPARC; 
Hanasoge 2010) code uses techniques developed by Hanasoge et al. (2006) and 



Hanasoge, Duvall, and Couvidat (2007) to simulate helioseismic wave propa- 



gation in the near-surface layers of the Sun. Waves are stochastically excited 
by introducing a forcing term in the vertical momentum equation; the forc- 
ing function is prescribed such that a solar-like power spectral distribution is 
obtained. The solution is temporally evolved using a second-order optimized 



Runge-Kutta integrator (Hu, 1996). The vertical derivative is resolved using 



sixth-order compact finite differences with fifth-order accurate boundary condi- 



tions ( Hurlburt and Rucklidge, 2000 ) . The derivatives in the horizontal directions 
are computed using FFTs (periodic horizontal boundaries) . The upper and lower 
boundaries are lined with damping sponges in order to enhance wave absorption. 
The SPARC code has been utilized to study wave propagation through model 



sunspots {e.g. 


Hanasoge 


2008 


Moradi, Hanasoge, and Gaily 2009 Moradi and 


Hanasoge 


2010 


1, as well as solar convection ( 


Hanasoge, Duvall, and DeRosa, 



2010a) 



A recent bit of progress with regards the choice of boundary conditions has 
been the development of a stable, unsplit PML formulation for the stratified 
linearized ideal MHD equations. Some related formulations have been developed 



by Parchevsky and Kosovichev ( 2007 ) and Khomenko and Collados ( 2006 ) . How- 
ever, instabilities caused by waves at grazing incidence to the boundary prevent 
long time integrations. By extending the technique of Convolutional Perfectly 



Matched Layers (C-PML; e.g. Roden and Gedney 20001, Hanasoge, Komatitsch, 



and Gizon (2010b) have succeeded in devising a stable C-PML formulation. 



6.3.4. The SAC Code 



In Shelyag et al. (2009), the propagation and dispersion of acoustic waves in a 



solar-like 2D sub-photosphere with localized, non-uniform magnetic field con- 
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centrations was investigated using the Sheffield Advanced Code (SAC) devel- 



oped by Shelyag, Fedun, and Erdelyi (2008). The numerical code is based on 



a modified version of the Versatile Advection Code (VAC; |Toth, Keppens, and 



Botchev 1998) , and employs artificial diffusivity and resistivity in order to stabi- 



lize the numerical solutions and relies on variable separation to background and 
perturbed components to treat gravitationally stratified plasma. The complete 
MHD equations are solved using a fourth-order central difference scheme for the 
spatial derivatives, and are advanced in time by implementing a fourth-order 
Runge-Kutta numerical method. 



The standard Model S ( Christensen-Dalsgaard et ai, 1996) atmosphere is 



employed as an initial background model, modified so that Fi is kept constant 
in such a way as to have the adiabatic sound speed profile closely match the 
sound speed in Model S. Three different self-similar, non-potential, magnetic- 
field configurations are employed for the simulations. As the curvature of the 
magnetic field changes the temperature stratification in the domain, the mag- 
netic configurations chosen differ by their field strengths and inclination at the 
visible solar surface. The acoustic source generates a temporarily localized wave 
packet with the duration of about 600 seconds, which has a peak frequency 
of approximately 3.33 mHz. The amplitude of the source is of the order of a 
few centimeters per second. Such a low amplitude makes sure that convective 
processes will not be initiated in the otherwise convectively unstable equilib- 
rium, and that the perturbations are kept linear, i.e., they do not change the 
background strongly. 

Three weak strongly-curved magnetic field with = 120 G at the 

surface, a strong weakly-curved magnetic field with Bz — 3.5 kG, and a strongly- 
curved strong magnetic field with Bz — 3.5 kG, were analyzed by means of 
time-distance helioseismology. The travel time dependencies show that for the 
first bounce the main part of effect of magnetic field on the acoustic wave is 
due to the change of the temperature structure in the sunspot. Nevertheless, the 
wave mode conversion from purely acoustic to the slow magneto-acoustic wave 
motion, characterized by an energy leak downwards, is also observed in the cases 
with strong surface magnetic field. 

6.4. Eikonal Methods 

While numerical simulations of wave propagation have significantly aided our 
level of understanding of helioseismology, further guidance is still needed in 
setting up the correct numerical experiments and understanding wave propaga- 



tion in magnetized plasmas. MHD ray theory (Weinberg, 1962) has traditionally 
provided a very useful conceptual framework in which to understand wave prop- 
agation, even though this is questionable at the surface where the pressure and 
density scales vary rapidly, since the assumption of slowly varying coefficients 
may not be justified. 

Regardless of these shortcomings, however, ray theory has been used in he- 
lioseismology for some time, being one of the several methods that have been 
applied to asymptotic inversions of helioseismic frequency measurements in the 



past {e.g. Cough 1984). In general, it has performed well beyond its formal 
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domain of applicability, a prime example being the agreement between the wave 
mechanical analysis of Cally ( 2005 ), the ray theory modeling of Cally ( 2006 ) and 
the recent results of Hansen and Cally ( 2009| ), who find very good agreement 
between generalized ray theory and previously published exact solutions ( |Cally,| 



2001 Cally, 2009b). 



Moradi and Cally (2008) recently combined eikonal methods and observa- 



tional data by constructing a 3D sunspot model based on observed surface 
magnetic- field profiles to propagate magneto-acoustic rays across the sunspot 
model for a range of depths to reproduce a skip-distance geometry similar to 
center-to-annulus cross-correlations used in time-distance helioseismology. In 
another recently completed work, Moradi, Hanasoge, and Cally (2009) com- 
pared the results of forward modeling via MHD ray theory and a 3D ideal-MHD 
solver, concluding that the simulated travel-time shifts were strongly determined 
by MHD physics and confirmed their strong dependence on the frequency and 
phase-speed filter parameters used. Khomenko et al. (2009) have similarly ap- 
plied MHD ray theory to validate and analyse the results of their numerical 
simulations. In another series of recently published works, Cally (2009a I investi- 
gates the direct (magnetic) and indirect (thermal) effects of the magnetic field on 
vertically propagating waves, suggesting that, overall, travel-time perturbations 
in umbra appear to be predominantly thermal, while in penumbrae they are 
mostly magnetic. Cally| (2009b) provides strong evidence for significant phase 
jumps (or discontinuities) associated with fast magneto-acoustic rays that pen- 
etrate the a ~ c level in sunspots. This effect appears to be more pronounced in 
highly inclined field characteristic of penumbrae. 

Most of the above work has centered on the study of properties of individual 
rays propagating through magnetized atmosphere by solving the ray equations 
for a point source in non-magnetic solar model. In a recent study, [Shelyag et al.\ 
(2009) have shown the importance of considering the full family of rays corre- 
sponding to a particular problem's initial conditions, as these define important 
geometric properties of the excited wave-field such as the wave front, caustics, 
and phase surfaces. 



7. Update on the Analysis of AR 9787 

7.1. Travel Times Comparison: Time-Distance and Helioseismic Holography 

It has been thought by some that travel times computed from time-distance 
helioseismology (Duvall et al, 19971 should be very similar (or maybe identical) 
to those measured using helioseismic holography (Braun and Lindsey, 2000). 
However, a detailed comparison of measured travel times has not been done. 
We have done such a comparison for an 8-hour interval for AR 9787. For the 
time-distance case, correlations are calculated between the central point and the 
quadrants. For the holography, egression and ingression signals are constructed 
in a quadrant geometry and correlations are done with the central point. Normal 
phase-speed filters were used (Couvidat et al., 2005|). 



Cuts across the travel-time maps for the sunspot are shown in Figure [6] for 
two distances. Frames a) and b) show the distance 24.35 Mm and c) and d) 
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show the shortest distance 6.20 Mm. The agreement is excellent for the distance 
of 24.35 Mm. However, there are still disagreements for the shortest distances, 
which we hope to understand soon. 
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Figure 6. Comparison of travel times computeci from helioseismic holography (red) with 
time-distance helioseismology (blue). Cuts through the center of the sunspot in the east- west 
direction are shown, a) Inward-going times, distance of 24.35 Mm. b) Outward-going times, 
distance of 24.35 Mm. c) Inward-going times, distance of 6.20 Mm. d) Outward-going times, 
distance of 6.20 Mm. 



7.2. Frequency Dependence of One- Way Travel Times 

The dependence of mean and difference travel times on the central frequency of 
the wave packets measured over a sunspot region was shown and interpreted by 



Braun and Birch ( 2006 2008 ) , as an indication of perturbations largely confined 
to a region not deeper than a few Mm. The largest frequency variations in travel 
times were seen for the small travel distances A, which are of the size (diameter) 
of the spot or smaller. In this section, we conduct a similar study for the sunspot 
in AR 9787. This type of study is important as, in principle, it should help us 
to constrain models of the subsurface structure of sunspots. 

Apart from a phase speed filter centered around a phase speed Vph~ 54 km 
s^^ with a width (FWHM) of about 40% of Vph, we also use Gaussian frequency 
filters, of 1-mHz width, centered at every 0.25 mHz interval between 2.5 and 5.0 
niHz to study the frequency dependencies of the travel times. In Figure[7]we show 
the umbral averages of perturbations associated with in-going and out-going 
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Figure 7. Frequency dependence of surface-focus travel times (shown as one-way and mean 
umbral averages of perturbations) for A = 50 Mm, for the sunspot in AR 9787. For clarity, 
only error bars for mean travel times are shown. The error bars represent standard deviations 
of travel time perturbations over the umbral pixels. 



wave travel times (measured using center-to-anmilus surface focus geometry), 
(Jrin and ^Tout, as well as mean travel-time perturbations, 5r„iean, against fre- 
quency for A = 50 Mm. There appears to be a strong frequency dependence 
associated with both Srin and Jrout measured across the umbra of the sunspot 
in AR 9787. This result appears to be in good qualitative agreement with the 



measurements of Braun and Birch (2008), who look at frequency dependence of 
the travel-time difference, albeit for another sunspoi[^ 

The underlying reason for such a strong frequency dependence of travel times 
at large distances is not entirely clear however. The interaction of helioseismic 
waves with sunspots can be influenced by a number of effects: the vertical ex- 



tent of the sunspot and the Wilson depression, p-mode absorption (e.g.. 


Braun, 


Duvall, and Labonte 


1987 


Braun 1995 


Cally, Crouch, and Braun 


20031, the 


contribution from subsurface flows, radiative transfer effects (Rajaguru et al, 



2006) and the impact of power reduction and source suppression in sunspots 
([Rajaguru et al, 2007 Hanasoge et al., 



2008 



Chou et al, 2009[ ). More detailed 
analyses and modeling of the sunspot in AR 9787 (see e.g. Section 7.8 1 will 



be required to isolate these effects and identify the cause(s) of the observed 
frequency dependence of travel times. 



Braun and Birch 



^In particular, see top right panel of Figure 5 in 
travel-time difference results for Filter I (distances of 48 - 55 Mm 



( 2008 I , which contains the 
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7.3. Effects of Filtering on Travel Times 

In this section we show a simple toy model for the sensitivity of time-distance 
travel times to changes in the power spectrum of solar oscillations. We will use 
this simple model to develop a qualitative understanding of the ridge and inter- 



ridge measurements (Thompson and Zharkov, 2008) of AR 9787 (Gizon et al, 
2009aj ). 



7.3.1. A Model Power Spectrum 

Before studying the impact of data analysis filters on time-distance measure- 
ments, it is necessary to have a model for the power spectrum of solar oscillations. 
Here we describe a simple model obtained from fitting the azimuthal average 
(over the angle of the wavevector k) of a 24 hr power spectrum from full-disk 
MDI data from a quiet Sun region. We then carry out a least squares fit to the 
azimuthally averaged power spectrum with a model spectrum of the form, 

P{k,u;)^\G{k,iu)\'+B{k,u;) , (14) 

where k — ||k|| is the horizontal wavenumber, uj is the temporal frequency, B is 
the background power, and the function G is given by 

G{k,.) = g «"(^)|[^_^„(fc)]^„(fc) - i^7TZ4M<ik)\ • ^'^^ 



This form of G is motivated by Equations (47) and (48) of Birch, Kosovichev, and 



Duvall ( 2004 \ in the case where the imaginary parts of the mode eigenf unctions 
are small. The free parameters in the fit are the complex mode frequencies a;„(fc), 
the mode amplitudes a„(fc), and the parameters describing a background power 
spectrum B that is linear in uj at each value of k. Here n is the radial order 
and rtmax is the maximum radial order used in the normal mode summation. In 
the examples shown here, Timax — 2. Throughout this toy model we will work in 
plane-parallel geometry; this is appropriate as we will be considering distances 
and wavelengths that are small compared to the solar radius. 

7.3.2. Sensitivity of Travel- Time Shifts to Change in Mode Frequencies 

For problems where the statistics of the wavefield are horizontally translation 
invariant and the expected (limit) power spectrum is azimuthally symmetric, the 
expectation value of the time-distance cross-covariance C(A, t) can be obtained 
from the limit power spectrum P(fc, a;) as 

/oo poo 
duj / kdk Ja{kA)T'^{k,uj)P{k,u)e-''^* (16) 
-oo Jo 

where Jq is the zero-order Bessel function, and J- is the data analysis filter 
(see e.g. Gizon and Birch 2002). Using Equation (16 1, we can compute the 
cross-covariance that we would expect for any model of the power spectrum. 



SOLA: SDLA1202R2_Moradi.tex; 24 August 2010; 1:25; p. 43 



H. Moradi et al. 




-50 50 

6(ii)/2jt (fxHz) 

Figure 8. Travel-time siiifts for the case of tlie ridge filter (thin lines) and inter-ridge filter 
(heavy lines) as functions of the change Slj in the real part of the resonance frequency of 
n = 1 mode. The solid lines show travel-time shifts obtained from the one-parameter method 
of Gizon and Birch (2002) and the dashed line show those from the linear definition of Gizon 
and Birch (2004). In both cases, the fitting window is chosen to be twenty minutes wide and 
centered at 28 minutes time lag. The travel distance in all cases is A = 16 Mm. Notice that 
for the case of the ridge filter, the two definitions of travel-time shift give very similar results. 
For the case of the inter-ridge filter, the difference between the two definitions increases with 
the amplitude of the frequency shift. 



Travel-time shifts can be obtained from the cross-covariance function using a 
wide variety of techniques {e.g. Duvall et al. 1997 Gizon and Birch 2002 Gizon 



and Birch 2004). For the examples shown here we use the linear definition of 



Gizon and Birch (2004) and also the one-parameter fitting method from Gizon 



and Birch (2002). 



Thompson and Zharkov (2008) showed that in sunspots, travel-time shifts 



measured with ridge filters (filters which isolate power along ridges) and inter- 
ridge filters (filters which isolate the part of the k — u! between the ridges) give 
travel-time shifts (relative to quiet Sun) of opposite sign. They found that the 
ridge filters yielded decreased phase travel times. This would seem to imply the 
waves in sunspots have increased phase speed. The same behavior was observed 



in the case of AR 9787 by Gizon et al. (2009a 2009b). 



As a very highly simplified toy model of this situation, we investigate the 
travel-time shifts caused by increasing the resonant frequencies in the power 
spectrum. The reason for choosing increased mode frequencies is to mimic the 
situation of increased phase speeds {i.e., increased mode frequency at fixed 
k). The procedure is as follows: i) use Equation (16 1 to compute a reference 
cross-correlation from the model power spectrum described above, n) compute 
a perturbed power spectrum by increasing the model resonance frequencies by 
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Figure 9. Slices through the filtered and unfiltered power spectra for the ridge filter (left) 
and inter-ridge filter (right). For the ridge filter case the slice is at fc = 1.1 radMm"^ and for 
the inter-ridge case the slice is at fc = 0.94 radMm"^ (in both cases the slice is near the k 
of maximum power). In both panels, the red (black) dotted lines show the unfiltered power 
spectra with (or without) the frequency shift. The green lines show the filter. The solid red 
(black) lines show the filtered power spectra with (or without) the frequency shift. In the ridge 
filter case, the consequence of increasing the mode frequency is that the mean power moves 
to higher frequency (thus higher phase speed). In the inter-ridge filter case the mean power 
moves to lower frequency (lower phase speed). 



a constant value (5w, Hi) compute a perturbed cross-covariance from this per- 
turbed power spectrum, and iv) fit tlie perturbed and reference cross-covariances 
to obtained the shift in the travel time caused by the change in resonance 
frequencies. 

Figure [8] shows the results of this procedure for two different choices of the 
filter T: the pi ridge filter and the inter-ridge filter isolating the region between 
the / and pi of Thompson and Zharkov (2008). In both cases a 1 mHz bandpass 
filter centered at 3.5 mHz has also been applied (Thompson and Zharkov, 2008). 
In the results shown here we have only considered the impact of changes in the 
frequency of the pi ridge. For the case of the ridge filter, increases in the mode 
frequencies (i.e., increased phase-speeds for all waves) yield decreased phase 
times. This is the expected result. For the case of the inter-ridge filter, however, 
we find that increases in wave speeds yield increased travel times. This is an 
unintuitive result and in this toy model is due to interaction of the inter-ridge 
filter with the power between the ridges. 

Figure |9] shows slices through power spectra for the two cases. The main 
effect of the perturbations to the mode frequencies is to move the ridges in the 
unfiltered power spectra. For the inter-ridge filtered case, the effect is more subtle 
as this filter isolates the part of the spectrum that lies between the ridges. The 
net result of moving the pi ridge to higher frequency and applying the inter-ridge 
filter is to move power to lower frequency (hence lower phase speed) at fixed k. 
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7.4. Ring-Diagram Structure Inversions 



The mode frequencies from ring-diagram analysis can be used to determine 
the thermodynamic structure of the gas under the tracked region. The mode 
parameters themselves are obtained by fitting a function P„(w, k^, ky) to the 3D 
power spectrum. 




5 10 15 20 

Depth (Mm) 



Figure 10. The effects of inversion parameter choices on inversion results. Three of the four 
inversion parameters are varied in these inversions - these inversions were performed with a 
fixed target width oiW = 0.0055. The top panel shows inversions with a surface term removed, 
and the bottom panels show inversions without a surface term. In each case, three different 
inversions are shown with different values of the error suppression term fi and the cross term 
suppression term /3. 



In this work, we use a fourteen parameter function defined in Basu, Antia, 



and Bogart (2004). The fit ridges can be interpolated to wavenumber k, and the 



resulting frequencies a;„(fc) can be treated as normal modes using the techniques 
of global mode analysis. Typically, inversions for solar structure are performed by 
linearizing the stellar oscillation equations around a reference model. Then, the 
differences in the frequency oJi of the ith mode between the data and reference 
model, with i representing the pair (n, k), are related to the changes in structure 



SOLA: SDLA1202R2_Moradi.tex; 24 August 2010; 1:25; p. 46 



Modeling the Subsurface Structure of Sunspots 



by the following equation: 



LU,, 



Kl.{z&^) dz + 



p Qi 



(17) 



The surface term, fgu^f (wi), is a smoothly varying function of frequency which 
accounts for non-adiabatic effects confined to the surface layers of the Sun, and is 
normalized by the mode inertia Qi. The observational errors are given by e^. The 
inversion kernels, K^^^iz) and ifp(z), are known functions of the reference model, 
and give the sensitivity of a mode to changes at a given depth. In ring-diagram 
analysis, however, we generally invert relative to mode frequencies measured in 
a nearby inactive region of the Sun. 




10 

Depth (Mm) 



20 



Figure 11. Ring diagram inversions for the structure beneath AR 9787. Inversions are per- 
formed for both sound speed squared (top panel) and adiabatic index Fi (middle panel). 
The inversions are performed relative to a nearby quiet region - the sense of the inversions 
is active minus quiet. The results for and Fi are used to infer a profile for the density p 
(bottom panel). 

The inversion technique used is called Subtractive Optimally Localized Av- 



erages (SOLA), a technique pioneered in terrestrial seismology (Backus and 



Gilbert, 1967). The SOLA method takes a specified target averaging kernel 
T{zo, z) which is localized around a target height zq, and minimizes the difference 
between that target and the actual averaging kernel, along with the contributions 
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from the cross-term kernel and from the errors (see e.g. Rabeho-Soares, Basu, 
and Christensen-Dalsgaard 1999 ) . There are two trade-off parameters, /it, which 



acts to suppress errors, and /3, which suppresses the cross-term kernel. There 
are also four free parameters to be chosen in this inversion method: the cross- 
term suppression term /3, the error suppression term fi, the characteristic width 
of the target kernel W, and A, the number of polynomials used to expand 
-Fsurf- Structure inversions using ring-diagram modes are quite unstable, and 
very sensitive to the choice of these parameters. Figure [TO] shows sound speed 
inversions with a variety of different choices of inversion parameters. The most 
dramatic choice to be made is whether or not to include a surface term in the 
inversions, and inversions are shown both with and without a surface term. A 
good inversion should have a well- localized averaging kernel )C{z(),z) around a 
target height zg, as well as reasonably small cross-term kernels. 

In Gizon et al. ( 2009a 2009b I, we presented an inversion of the region AR 9787 
for sound speed. In Figure 11 we show an inversion of the same region for sound 
speed and for adiabatic index Fi. The inverted quantities are fairly consistent 



with the results from other active regions {e.g. Basu, Antia, and Bogart 2004 



Bogart et al. 2008 1 . There is a depression in sound speed between approximately 



3 Mm and 10 Mm depth. Below that depth, the perturbation becomes positive. 

The adiabatic index is depressed below 3 Mm, with an enhancement below 
12 Mm. If and Fi are determined, the other thermodynamic quantities are 
also, in principle, known. In Figure [TT] we also show an inferred density profile. 
The density between approximately 5 Mm and 11 Mm appears to be slightly 
depressed in the active region compared to the quiet Sun, while the region below 
that has a lower density than the quiet Sun. There is a depression in density at 
the shallow edge of the inversion, which might be a sign of a Wilson depression 
of the optical surface. It should be noted here, however, that results in the 
shallowest layers of the Sun are quite uncertain, due both to lack of resolution 
in the inversions and to uncertainties in the physics. 

7.5. Moat Flow: Hankel Analysis 

Fourier-Hankel decomposition is a useful analysis tool to study p modes in 
the vicinity of sunspots. In particular, the method has been used to study the 



absorption of p-mode power by active regions (Braun, Duvall, and Labonte, 
19881 iBogdan et al., 19931 iBraun, 19951 IChen, Chou, and TON Team, 19961). In 



this procedure, the wave signal of the p modes is decomposed into inward and 
outward propagating modes in an annular region surrounding a sunspot. The 
annulus is chosen to be small enough to allow describing the spatial part of the 
solar oscillations with Hankel functions. In the seismic study of AR 9787 un- 



dertaken by Gizon et al. (2009a 2009b ) , the power absorption of this particular 
sunspot was demonstrated by Fourier-Hankel decomposition for the m = p 
modes. 

Fourier-Hankel decomposition also turns out to be useful for investigating 
the horizontal outflow associated with active regions. An outflow directed hori- 
zontally (the moat flow) would result in Doppler shifts of the p modes traveling 
into and out of the spot. A radial flow therefore leaves its signature in the power 
spectra by shifting the p-mode ridges. 
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Figure 12. Power spectra for two inward (solid line) and outward (dotted line) propagating 
p modes. Left: £ = 288, n = 3; right: £ = 452, n = 2. The peaks of the outward propagating 
modes are shifted to higher frequencies. 



Here we extend the Fourier-Hankel analysis of AR 9787 using the procedure 



described by Braun ( 1995 1. For this analysis, the annular region around AR 9787 
is defined by an inner radius of 30 Mm and an outer radius of 137 Mm measured 
from the spot center. Power spectra for inward and outward propagating waves 
were obtained for modes with azimuthal order m = —5, —4, . . . , 5. Figure [12] 
displays the resulting spectra for two modes. It is immediately apparent by 
visual examination of the spectra that the outgoing p-mode power is shifted 
to higher temporal frequencies relative to the ingoing p modes. We note that 
the power spectra displayed were normalized in order to correct for the power 
absorption due to the presence of the sunspot. The observed frequency shift is 
of the order of 10 /iHz. 

7.6. Ring-Diagram Analysis of Flows 

7.6.1. Large Scale Flows Around Active Regions 

Dense-pack ring-diagram analyses of MDI and GONG-I--I- data have shown that 
subsurface flows associated with active regions have the following characteristics: 

• Active regions are surrounded by extended inflows and outflows, depending 
on depth. 

• Subsurface flows of active regions are twisted. 

• Active regions rotate faster than quiet regions. 

• Flux emergence correlates with the flow divergence signal. 

The horizontal components of solar subsurface flows are determined over a range 
of depths from the surface to about 16 Mm using the dense-pack ring-diagram 



analysis (Haber et al., 2002). Daily flow maps are calculated for 189 dense- 
pack regions of 15° diameter with centers spaced by 7.5° in latitude and central 
meridian distance. The dense-pack technique thus measures flows on horizontal 
scales comparable to the size of active regions. 
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Locations of strong active regions show, on average, extended divergent hori- 
zontal flows at depths greater than about 10 - 12 Mm and convergent horizontal 
flows closer to the surface (both with an amplitude of about 50 m s~^). AR 9787 



shows this pattern as well (Figure 13 ). Using a mass conservation constraint, ring 
analysis infers vertical flows associated with the horizontal flows. The conver- 
gent and divergent flow pattern, is one of the most consistent characteristics of 
subsurface flows associated with active regions. It has been studied with several 



2004 



helioseismic techniques (Gizon, Duvall, and Larsen, 2001 Braun, Birch, and 



Lindsey, 2004 Haber et al., 2002 Haber et al., 2004) Zhao and Kosovichev, 



Komm et al., 2005 ). These extended flows, however, should not be confused 
with the much more localized moat flow (divergent, 250 m s~^, see next section 



and Gizon et al. 2009a, 2009b), which cannot be resolved with 15° aperture ring 
analysis. 

Active regions are locations of flows with larger vorticity values; they are 
locations of strong vertical gradients of the horizontal flows. Strong active regions 



are identiflable by a dipolar pattern in zonal and meridional vorticity (Mason 



et al, 20061 and AR 9787 appears to be no exception (Figure 13). The presence 
of active regions is barely noticeable in vertical vorticity maps of this spatial 
resolution. However, active regions are, on average, characterized by cyclonic 
vorticity (counter-clockwise in the northern hemisphere), which might be due to 



2003 Zhao and Kosovichev, 2004 1. 



the Coriolis force acting on the flows (Spruit, 2003). This agrees with observa- 



tions with higher spatial resolution ( Duvall and Gizon, 2000 Gizon and Duvall, 



From direct surface measurements, it is well known that active regions rotate 



faster than quiet ones. This has been conflrmed with helioscismology (Gizon, 



2004 Zhao and Kosovichev, 2004). A ring-diagram analysis of about six years of 
GONG-I— I- data shows that the average zonal flow of active regions is about 4 m 
s~^ larger than that of quiet regions from the surface to a depth of 16 Mm (Komm 



et al., 2009b). The difference is about one order of magnitude smaller than that 
derived from surface measurements of active and quiet regions, which is most 
likely a consequence of the rather large size of the dense-pack patches and the 
resulting averaging over many different types of magnetic features. Results from 
acoustic holography and time-distance analysis with higher horizontal resolution 



support this interpretation ( 


Braun, Birch, and Lindsey, 2004 


Gizon, 2004 


Zhao, 


Kosovichev, and Duvall, 2004 


I- 



A survey of 788 active regions observed with GONG-I--I- makes it possible to 
determine a signature of emerging magnetic flux in subsurface flows associated 



with active regions (Komm, Howe, and Hill, 2009a). At depths greater than 



about 10 Mm, upflows become stronger with time when new flux emerges. 
At layers shallower than about 4 Mm, the flows might start to change from 
downflows to upflows, when flux emerges, and then back to downflows after the 
active regions are established. The flow response to emerging flux agrees with 



numerical simulations of emerging flux tubes (Fan, 2001 Schiissler and Rempel, 



2005 ) where upflows indicate the beginning of flux emergence and surface cooling 



due to adiabatic expansion leads to downflows along the emerged loops. 
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Figure 13. Top: The unsigned magnetic flux (solid line) at 127.5° — 135.0° longitude as a 
function of latitude and binned over 15° (dotted line). Bottom: The kinetic helicity density at 
127.5° - 135.0° longitude as a function of latitude and depth. The kinetic helicity density is the 
scalar product of the velocity and vorticity vector. The arrows represent the meridional and 
vertical velocity components with the vertical one increased by a factor of ten for visibility. 
Average zonal and meridional flows have been subtracted. Active region 9787 is noticeable as 
the location of strong helicity values of opposite sign which coincides with the peak of the 
unsigned magnetic flux. The region shows upflows at depths greater than about 10 Mm and 
downflows at shallower depths at grid points of -15° to -7.5° latitude. 



7.6.2. Sunspot Flows from Small Rings 

Subsurface horizontal flows determined by high-resolution ring analysis (HRRA) 
for the four days when AR 9787 was closest to the center of the disk, show 
characteristic outflows from the lone sunspot corresponding to moat flows (see 



Figure 14 1. 

The analysis was carried out on 2°-diameter tiles whose centers were spaced 
1° apart. Each flow arrow represents a spatial average over an entire tile and is 
itself the result of averaging the flows determined from all the fitted / modes of 
the power spectrum for that tile since there are not enough modes to perform a 
true inversion. This means that the flows are characteristic of the gas from the 
surface down to a depth of to 2 Mm where the / modes reside. The four-day 



panel (Figure 14 ) shows the evolution of large-scale zones of divergence around 
the sunspot corresponding to large supergranules as well as a seeming twist of 
the flows coming from the sunspot on 25 January. The flows are over-plotted on 
averaged MDI magnetograms for the given day where the green and red colors 
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Figure 14. Subsurface flows using high-resolution ring analysis of full-disk MDI Dynamics 
Doppler data from 22-25 January 2002. MDI magnetograms are underlaid with red and 
green specifying opposite polarities of the magnetic field. The flows were determined from 
/-mode data and are thus representative of flows within the top 2 Mm of the convection zone. 
The outflows from the sunspot (shown in green) correspond to moat flows, while the cells of 
divergence seen in the vicinity of the active region correspond to larger supergranular scales 
of about 50-60 Mm. The magnetic network often appears at the edges of these larger cells. 



show magnetic fields of opposite polarity (the green circle corresponding to the 
sunspot). 



Similar subsurface flow signatures were obtained by Gizon et al. (2009a 



2009b), who used / to p4 ridge-filtered time-distance travel times to produce 



linear inversions for flows (see Jackiewicz, Gizon, and Birch 2008) around the 
sunspot in AR 9787. Figure[T5]shows a summary plot of the subsurface horizontal 
flows around the sunspot, which appear very much consistent with the flows 
detected by HRAA, and direct observations of the moat flow (see Section 2.2.5 1. 
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Figure 15. Map of horizontal flows (arrows) at a depth of 1 Mm around the sunspot in AR 
9787 using one day of MDI full-disk data and time— distance helioseismology on January 24 
(Gizon et al., 2009). The spatial resolution is determined by the width of the averaging kernel 
(FWHM 7 Mm) shown in the top-left corner. The longest arrow corresponds to a flow of 450 m 
s~^. The surface line-of-sight magnetic field is displayed in red and blue shades (saturated at 
±350 G). This inversion uses the modes / and pi through p4; it is discussed by Jackiewicz, 
Gizon, and Birch (2008). 



7.7. Acoustic Halos 



7.7.1. Observations 



The acoustic halo is an observed enhancement of high frequency (i.e., above the 
acoustic cut-off frequency at approximately 5.3 mHz) acoustic power surrounding 
regions of strong magnetic field. Hindman and Brown ( 1998 ) demonstrate that 



this enhanced power at high frequencies tends to be prominent in intermediate 
magnetic-field strengths of 50-250 G and appear to be absent in the equivalent 
continuum intensity observations. In Gizon et al. (2009a, 2009b) it was demon- 



strated that AR 9787 has an extensive plage region, with the excess power at 
above 5.5 mHz being spatially correlated with these plage regions. It was also 
noticed that the strongest regions of acoustic halo occur between regions of 
opposite polarity. 

For our analysis, we calculate the vector magnetic field at the surface from the 
MDI magnetograms using a potential magnetic-field extrapolation. From this we 
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define the inclination of the magnetic field from vertical [7]. Figure 16 shows 
the average normalized power against total magnetic-field strength for 7 < 45° 
(black) and 7 > 45° (blue) for 1 mHz frequency band passes centered at 3, 4, 5 
,and 6 mHz as indicated. This highlights the importance of the field inclination 
for the acoustic halo. 



In another recent analysis, Schunker and Braun (2010) further demonstrate 



that the acoustic halo is specifically confined to intermediate field strengths 
(between 100 G and 400 G) with 60° < 7 < 120° from vertical. The power 
confined to this range of inclination has a peak which increases in frequency 
from 5.5 mHz to 6.5 mHz as the magnetic-field strength increases from 100 G to 
400 G. In addition, they also find that the radial order ridges of the azimuthally 
summed power spectrum are shifted to higher wavenumbers than the quiet- 
Sun at constant frequency in the halo regions. The reason for this is, as yet, 
unexplained. 

7.7.2. Theories 

A growing number of theories attempt to explain the phenomenon of wave veloc- 
ity enhancements in the vicinity of active regions. The hypothesis of an altered 
wave excitation mechanism due to the presence of magnetic fields was proposed 



by Brown et al. (19921 and Braun et al. (1992), and revisited more recently by 
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Jacoutot et al. (20081. Of course, it is also equally likely that the substantive 



differences in the wave propagation physics in active regions could also be leading 



to these enhancements, a possibility noted by a number of authors {e.g. Hindman 



and Brown 1998 Donea, Braun, and Lindsey 1999 Moretti et al. 20071 



From analyses of simulations of waves propagation through a model sunspot, 
Hanasoge' f2008^ demonstrated that it was possible to obtain these halos without 



the requirement of enhanced sources. A wave scattering theory based on pure 



hydrodynamics was also put forth by Kuridze et al. (2008 2009); their thesis 



was that bi-polar canopies create a secondary trapping region for the outward 



propagating high-frequency waves. Alternately, Khomenko and Collados ( 2009 ) 
have suggested that upward propagating high frequency waves undergo conver- 
sion to fast modes, refract off the large Alfven speed gradient in the atmosphere, 
and return to the photosphere. This results in the same wave being observed 
twice, once while going upwards and the second time, on its way back down. 



thereby causing an enhancement. Finally, Hanasoge ( 2009 ) has proposed that the 



presence of a large number of wave scatterers in the vicinity of active regions may 
cause preferential scattering into low mode-mass waves, i.e. waves whose energy 
is focused in the near-surface layers. Relative to the original (quiet Sun) set of 
modes, this scattered configuration contains stronger surface velocity signatures 
{i.e. lower net mode mass), thereby leading to the halos. Observationally, it can 
be confirmed that the velocity enhancements are primarily present at high-/, 
in the region of the power spectrum where the modes possess the lowest mode 
masses. 

In reality, it may well be that there are multiple mechanisms at play, i.e., 
altered wave excitation, reflection/refraction of waves by magnetic fields, and 
preferential scattering into waves with low mode masses. Further studies are 
required before any solid conclusions can be drawn. 



7.8. Comparison of Observed Cross-Covariances with Simulations 



Cameron, Gizon, and Duvall (2008) and Gizon et al. (2009a 2009b) showed 



that the cross-covariance of the MDI Doppler velocity is a very useful quantity 
to study the scattering of solar waves by the sunspot in AR 9787. The cross- 
covariance between two points is closely related to the Green's function between 



these two points (see, e.g. Gizon, Birch, and Spruit 2010 and references therein). 



By extension, the cross-covariance between the signal averaged over a line and 
any other point can be used to study the interaction of plane wave packets with 
the sunspot. 

The left panels in Figure [17] show such cross-covariances at particular values 
of the correlation time-lag, after the wave packets have traversed the sunspot. 



While Gizon et al. (2009a 2009b) had studied only /-mode wave packets, here 



we also show the observed cross-covariances for pi and p2 wave packets. The wave 
packets are constructed using standard ridge filters (wave power peaks near 3 
mHz). The observed cross-covariances show that, in all cases, solar waves speed 
up through the sunspot and that their amplitudes are reduced relative to the 
quiet Sun. The stochastic noise in the observations is reduced by exploiting all 
available symmetries, including the near-cylindrical symmetry of the sunspot. 
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Figure 17. Left panels: Observed MDI cross-covariance functions between the Doppler ve- 
locity averaged over the line x = -40Mm and the Doppler velocity at any other point {x,y) in 
the map. The sunspot AR 9787 is located at the origin and the circles show the boundaries of 
the umbra and penumbra. The cross-covariance is averaged over seven days. The top left panel 
shows the cross-covariance for the /-mode ridge at time-lag t = 185 minutes, the middle left 
panel for the pi-mode ridge at time-lag t = 165 minutes, and bottom panel for the p2-mode 
ridge at time-lag t = 145 minutes. The right panels show SLiM simulations of the propagation 
of /, pi, and p2 wave packets in the -|-x-direction through the semi-empirical sunspot model 
of Cameron et al. (2010). The vertical velocity in the simulation is directly comparable to the 
observed cross-covariances. 
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The cross-covariances also provide very strong observational constraints to 
test the validity of sunspot models. Here we have used the SLiM code to test the 
semi-empirical model of AR 9787 by Cameron et al. (2010) (see Section 4.3.1 
and Figure 4) by propagating /, pi, and p2 wave packets through the sunspot 
model. The initial conditions of the simulations are chosen such that the vertical 
component of velocity in the simulation is directly comparable with the cross- 
covariance function. As can be seen in Figure [iTj the simulations reproduce the 
basic features of the observations. This indicates that the sunspot model has a 
seismic signature that is close to the real sunspot, which is very encouraging. 



8. Discussion and Perspective 

8.1. Sunspot Structure: A Critical Assessment of Existing Models 

Recent years have seen an increased interest in questions relating to the structure 
of sunspots, fuelled by a fortunate coincidence of several factors. After decades 
of more gradual improvements, classical (visible light) observations have quickly 
improved in quality over the past few years. With the Swedish 1-m Solar Tele- 
scope, spatial resolution achieved has made a jump (with a resolution of 0.1" 
achieved regularly in the blue), while time coverage has vastly improved through 
the long time sequences from the Hinode satellite. Secondly, helioseismic obser- 
vations have now approached the point where they can be turned into powerful 
diagnostics of spot structure, as discussed elsewhere in this text. Finally, and 
perhaps most dramatically, the realism achieved by 3D radiative MHD simula- 
tions has now opened the perspective of a definitive physical interpretation of 
the phenomenology of sunspot structure, including umbral dots, light bridges, 
penumbral filaments, the Evershed fiow, the moat flow, and their relations to 
each other. It is likely that understanding will soon replace the historically 
evolved patchwork of mutually inconsistent views and physically dubious ad- 
hoc models. The following gives a brief perspective on these developments. For 



(2010) 



more detail, see the reviews in Scharmer (2009) and Nordlund and Scharmer 



As discussed in Section[3] the structure of a sunspot as observed at the sur- 
face is kept together by forces deeper down. Assuming this anchoring at depth 
as given (but yet to be explained), one can seek a physical interpretation of 
the surface structure observed in umbra and penumbra. At the observed field 
strength, 1-3 kG, the magnetic field is strong enough to suppress convection. 
This explains why spots are dark but it does not explain the particular time 
dependent fine structure observed. 



In early observational work {e.g. Mamadazimov 1972) the mix of long dark 



and bright filaments has been interpreted as showing a magnetic field (dark, 
as in the umbra) on top of the normal photosphere (shining through in the 
bright filaments). This explained the general appearance of the penumbra and 
the nearly photospheric brightness of the bright structures, but required that the 
penumbral magnetic field only touches over the photosphere (a "thin penumbra" , 
with optical depth only of order unity). The field would have to be essentially 
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horizontal in the penumbra, field lines crossing the solar surface only in the 
umbra. This does not agree with the observed inclination of the penumbral 
field. In fact, most of the magnetic flux of a spot crosses the surface through 
the penumbra, not the umbra. The penumbral field, on average, is therefore 
not horizontal, and the observed field must continue to some depth below the 
surface. 

This has led to interpretations in terms of convection in a magnetic field 
extending to a substantial depth below the surface (a "thick penumbra"). An 
influential conceptual picture was the Danielson (1961) model of convective 



"rolls" : an overturning flow in a plane perpendicular to a horizontal magnetic 
field. This idea has led to a "magnetoconvection" view of the penumbra, which 
interprets the observed filamentary structures as turbulent fluctuations in a mean 



magnetic field inclined at some angle to the solar surface (see e.g. Tildesley and 



Weiss 2004 and references therein). It also plays a role as a conceptual picture 
underlying global models of sunspot structure, such as Jahn and Schmidt ( 1994 ). 



A key aspect of any realistic model of magnetic surface structure is the 
transition from the gas pressure dominated regime below the photosphere to 
the magnetically dominated atmosphere. In the penumbra this takes place just 
at the photospheric level. The transition takes place over a couple of pressure 
scale heights. This is no more than the horizontal resolution of most telescopes, 
and comparable to the smallest horizontal structures seen in the penumbra. 
With photospheric observations one is thus looking at the thin interface between 
two physically different regimes, rather than a slice through some quasi-uniform 
turbulent medium. The consequences of this fact have not been realized in most 
of the currently popular views. 

Observations of this interface show a mix of magnetic-field strengths and 
inclinations. This has led to the idea of inclusions ("flux tubes") embedded in 
a background field of different direction (Solanki and Montavon, 1993). This 



has become the dominant theme in interpretations of penumbral structure {e.g. 
Martinez Pilietl [2000l [Bellot Rubio et al.\ [2003bl [Bellot Rubio, Balthasar, and| 
CoUadosI |2004[ |Borrero| |2007[ ). 

In the atmosphere, however, where B'^/8tt ^ P, forces other than magnetic 
are small. The magnetic field in such a case is nearly force free. This provides 
strong restrictions on the physically realizable field configurations. These con- 
straints are actually violated in most proposed ideas about penumbral structure, 
especially the floating flux tubes in the embedding and uncombed type models. 



An exception is the proposal by Martens et al. { 1996 ) , who propose a structure 



of magnetic sheets of different directions separated by force-free currents. 

Force-free, but not current-free, configurations also cause problems however. 
The strength of the magnetic field decreases away from the spot. A force- free 
field has the property that a difference in field-line direction between neighboring 
magnetic surfaces increases with decreasing field strength along the field (this is 
related to the magnetic torque being conserved along field lines, see discussion 
in Spruit and Scharmer 2006). Force free models thus predict that differences in 



inclination should become more prominent with height above the photosphere 
and distance from the umbra. Observations in the chromosphere (the "super- 
penumbra") do not support this. Instead, magnetic-field line directions traced 
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by chromospheric indicators appear to become more uniform with height. This is 
not consistent with force-free models hke those of Martens et al. (1996). Instead, 
these observations indicate that the field in the atmosphere is much closer to a 
potential (current free) configuration. 

The "gappy penumbra" model (Spruit and Scharmer, 2006 Scharmer and 



Spruit, 2006), based on the cluster model of a sunspot, explains how the pattern 



of strong irregularities decreasing with height above the surface comes about 
naturally in a potential field, as a result of gaps between field lines created below 
the surface by overturning convectiorj^ At the same time this model addresses 
the long-standing heat flow problem: the question how the observed radiative 
energy flux is supplied to the penumbral surface. The low field strength in the 
gaps allows convection to supply the heat flux emitted by the surface unimpeded 
by magnetic forces. Convection in magnetic rolls, on the other hand, docs not 
solve this problem, since a roll with the diameter of a penumbral filament does 
not contain enough thermal energy to maintain the radiative fiux over the life 
of a filament. However, the gappy model has problems explaining the structure 
in the outer penumbra, where observations show nearly horizontal field (also 
confirmed in the numerical simulations of'Rempel et al. 2009b ) , and the presence 
of returning magnetic flux (see Section 2.2.1 ). 

Apart from the gappy model, a number of other physical models for the 
penumbral currently exist. The model of [Schlichenmaier, Jahn, and Schmidt] 
( |1998j ) (see also Schlichenmaier 2009 and references therein), which has the 



virtue of substantial physical detail, prescribes that magnetic fiux tubes rising 
from below the surface carry heat and as well as an outward flow, interpreted as 
the Evershed flow. It has problems reproducing the observed heat fiux from the 
penumbra, however, and the postulated rising fiux tubes have not been identified 
in the numerical simulations. In the numerical simulations, both the heat flux 
and the Evershed flow are found to be driven by overturning convection, not in 
the form of magnetic rolls or flux tubes, but in regions low field strength (see 



Nordlund and Scharmer 2010). 



Another proposal appeals to the process of "turbulent flux-pumping" ( Thomasj 



et al., 2002 Weiss et al., 2004). In this process asymmetry between upflows and 



downflows in stratified convection effectively transports horizontal magnetic- 
field lines downward. Penumbral structure is then ascribed to this process. The 
source of the penumbral structure is said to be located outside the penumbra 
(Weiss et al., 2004). Alternatively, it has also been implied to be located in 
the penumbra itself (see Tildesley and Weiss 2004). The turbulent pumping 
mechanism appears to offer a plausible mechanism for producing the returning 



By "overturning" we mean here the convectiv e pattern also ob served in realistic simulations 
of stellar surface convection (see the review by |Scharmer| |2009| for more detail). As opposed 
to traditional views approximating convection as closed cells or "rolls" , almost nothing of the 
descending part of the flow returns to the surfac e but is replace d by upflows from larger depths. 
As opposed to the field-aligned rolls of the Danielson' ('1961') kind, this pattern maintains a 
low field strength in the gaps through the process of convective expulsion. While these gaps 
may not be exactly field free (as none of the convection zone is), they are regions where the 
hydrodynamic forces dominate over magnetic forces, i.e., the field is "below equipartition with 
convection" . 
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magnetic flux in the outer penumbra, as well as the apparent hysteresis observed 
in the transition between a pore and a sunspot (Weiss et al., 2004). However, 
there are also objections to this model that are based on the fact that the action 
of convection on a magnetic field overlying it, as in the outer penumbra, would 
be expected to have the opposite effect of turbulent pumping. Overturning con- 
vection would instead expel a magnetic field that is imposed from above, through 



the well-known process of convective expulsion (Zel'dovich, 1956 Weiss, 1966). 



The process is seen in action in numerical simulations of horizontal magnetic 



fields overlying granulation by Steiner et al. ( 2008 2009 1 . 



As the quality of observations continues to improve, realistic 3D radiative 



numerical MHD simulations of sunspots ( Schiissler and Vogler, 2006 


Heinemann 


et a/., 20071 'Scharmer, Nordlund, and Heinemann, 2008| Rempel, Sc 


liissler, and 


Knolker, 2009a, .Rempel et al, 2009b 


) show a remarkable level of agreement with 



these observations. The agreement includes the properties of umbral dots, the 
inward propagating bright filaments, the dark cores overlying them, the varying 
aspect of penumbral structure with viewing angle, the varying field strengths 
and direction in penumbral filaments, the dependence of these on height, the 
moat flow and the Evershed flow. We expect that more advanced simulations 
in the near future will further improve our understanding of penumbrae, in 
particular with regards to their outermost parts. With this level of agreement 
with observations, there is little doubt that the simulations are reproducing the 
physics of sunspot structure as observed at the surface. 

Next to the treatment of the atmospheric magnetic flcld, the physics of ra- 
diation is of equal importance for realism in numerical simulations. Cooling by 
radiation at the surface determines the thermal structure of the penumbra and 
drives the observed flows. On the other hand, it also determines the detailed 
appearance of penumbral structure at the optical depth unity surface. Any 
physically meaningful comparison with observations thus requires inclusion of 
radiation physics at a fairly well developed level. The fact that the level of realism 
needed for a meaningful interpretation of solar surface structure, magnetic or 
quiet, has now been achieved should be considered the most significant advance 
in theoretical solar physics of the past few decades. 



.2. The Starspot Connection 



It is important to keep in mind that our Sun is not the only star that possesses 
spots. Moreover, it is more than likely that all late-type stars with convective 



envelopes exhibit spots, or "starspots" {e.g. see reviews by Strassmeier 2009 and 
Berdyugina 2005 1 . Spotted stars constitute roughly 90% of all stars in the Milky 



Way, basically all GKM stars and a large fraction of F and L and T dwarfs are 
spotted, representing a mass range from the brown dwarf limit up to the 2.4 M© 
of a G giant, and an age range from the pre-main sequence phase up to the 
asymptotic giant branch. Stars with planets can also be affected by magnetic 
processes and their magnetic environment may even affect close-in planets and 



back- react on to the star {e.g. 


Shkolnik, Walker, and Bohlender 


2003 


Catala 


et al. 


2007 


Kashyap, Drake, and Saar 


2008 


Lanza 


2008 


). Therefore, we must 



understand sunspots first, before we can understand starspots. 
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Figure 18. A starspot on the K subgiant II Pegasi. The surface temperature and the photo- 
spheric magnetic field is a Zeeman-Doppler observation while the interior shows a numerical 
simulations of a non-axisymmetric dynamo in a fully convective star. There, polarity in the 
western hemisphere mirrors the polarity in the eastern hemisphere. On the surface, red means a 
field strength of up to 800 G according to the color bar while blue means basically no detection. 
In the interior, red means the field is pointing outwards, blue means the field points inwards. 



Models of starspots do not exist yet but the mere size difference suggests that 



a scahng from a sunspot model is not appropriate (see e.g. Solanki and Unruh 
2004). While sunspots typically cover 10~^ to 10~^ of the solar surface and only 
during solar maximum reach about 10~^, the record holder among other stars 
is still the one big starspot seen on the KO giant XX Tri in December 1997 
(Strassmeier, 1999). Based on the star's Hipparcos distance of 197 pc, the spot's 
area covers approximately 12 x 20 i?Q, i.e. 22% of the star's hemisphere or 10 000 
times the area of the largest sunspot group. Clearly, its emergence, structure, 
and decay must be addressed on a global scale. Figure 18 shows a Zeeman- 
Doppler image of the active spotted star II Pegasi (Carroll et al., 2007) combined 
with a numerical simulation of a non-axisymmetric dynamo. Although just a 
demonstration example, it nonetheless highlights the expected global starspot- 
dynamo connection. 
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8.3. Conflicting Helioseismic Observations 



A prevalent, largely phenomenological, approach to modeling the subsurface 
structure of sunspots has been to treat the regions of magnetism as pertur- 
bations to the background wave speed. These types of models have been con- 
structed using a variety of local-helioseismic procedures. A comparison of struc- 
tural (wave-speed) inversions for AR 9787 using both ring-diagram analysis and 
time-distance helioseismology with phase-speed filters was presented by |Gizon| 



et al. { 2009a 2009b I , and is partially reproduced in Figure 19 
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Figure 19. Comparison of different helioseismic metfiods used to infer wave speed perturba- 
tions below AR 9787 (Sc^/c). The red and blue curves show the ring-diagram and phase-speed 
filtered time-distance results, respectively, from Gizon et al. (2009a, 2009b). The time-distance 
result is shown along the axis of the sunspot; the ring-diagram results has been scaled by a 
factor of 10. The black curve indicates the fast mode speed perturbation ((cf — c)/c) from the 
radiative MHD simulations of Rempel et al. (2009b), described in Section [4.4| it approaches 
the value 75 at 5; = Mm. The green line represents the on-axis wave-speed perturbations 
deduced from the phenomenological model of Fan, Braun, and Chou (1995), based on the 
Hankel phase-shift measurements of Braun (1995). 



As is evident, these two inversions give subsurface wave-speed profiles with 



opposite signs and different amplitudes. Gizon et al. (2009a 2009b) discussed a 



number of factors which could contribute to such a disagreement, including the 
observation that the ring-diagram inversions include a treatment of near-surface 
effects absent from the time-distance analyses, and the apparent sensitivity of 
the measurements to details of the analysis which are not included in the mod- 
els. A prominent example is the sensitivity of time-distance measurements to 
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parameters of the filters (Section 7.3 1. As suggested by Braun and Birch (2008), 
some evidence of a strong near-surface contribution to hehoseismic measurements 
in sunspots is shown in the forward modehng of Fourier-Hankel measurements 



performed by Fan, Braun, and Chou (1995). It is worthwhile, for comparative 



purposes, to therefore plot the results of the subsurface perturbation suggested 



by Fan, Braun, and Chou ( 1995 1 in Figure 19 An important caveat in this com- 



parison is that the Hankel phase-shift measurements ( |Braun, 1995 1 , from which 



this wave speed result was inferred, were of different sunspots. However, Braun 



and Birch ( 2008 ) show that these phase-shift measurements agree favorably with 



more recent ridge-filtered holography measurements of sunspots of similar size. 
The relative fast-wave speed perturbations of the sunspot model of |Rempel et al.\ 
(2009b) presented in Section 4.4 are also included in Figure 19 for reference: this 



model also displays positive wave-speed perturbations in the first 2 Mm below 
the surface. 



8.4. Emergence of a New Paradigm in Sunspot Seismology 

Keeping in mind all of the above caveats, it is worth noting that three out 
of four curves shown in Figure [19] are consistent with a strong, positive wave- 
speed perturbation extending about 2.5 Mm below the surface. Below this depth, 
the hehoseismic inversions show considerably stronger deep wave-speed pertur- 
bations (albeit, of opposite signs) than the other methods. Braun and Birch 
(2006) have argued that shallow wave-speed perturbations like those suggested 



by Fan, Braun, and Chou ( 1995 ) are required to explain the systematic frequency 



dependence in the mean travel times observed in sunspots. Of course, it is well 



to keep in mind the warnings of Gizon et at. ( 2009a 2009b I , namely the possible 



naivety of modeling potentially complicated effects of the magnetic field in terms 
of an equivalent sound-speed perturbation. 

Many of these magnetic effects are more suitably explored through the numer- 
ical methods discussed in Section|6] Perhaps the strongest argument in favour 
of a shallow, fast wave-speed model is provided by the linear simulations of 
MHD wave propagation (Section [7.8[ ). These forward numerical computations 
show that a simple semi-empirical sunspot model extending no deeper than 
2 Mm (Figure [s]) is capable of reproducing many features of the hehoseismic 
measurements, in particular the cross-covariance signatures of the sunspot in 
AR 9787. While additional linear simulations will be needed to confirm this 
claim, the realistic radiative simulations of sunspot-like structure by [Rempel] 
et al.\ (|2009b|) will provide the ultimate testbed to validate the forward and 



inverse methods of sunspot seismology. 

Finally we note that sunspot seismology will benefit greatly from improved ob- 
servations by the Solar Dynamics Observatory (SDO), scheduled to be launched 
shortly. The Hehoseismic and Magnetic Imager (HMI) instrument onboard SDO 
will not only provide higher-resolution Doppler images of the solar disk, but 
it will also provide full vector magnetograms. Reliable measurements of the 
photospheric magnetic field will be key to constraining the near-surface layers of 
sunspot models, as detailed models of the surface layers are a necessity in order 
to probe the deeper structure of sunspots. 
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